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In  this  report  we  focus  on  the  roles  of  models  of  flexible  structures  in 
the  design  and  evaluation  of  control  laws  for  the  damping  of  vibrational 
motions  in  those  structures. 

In  the  first  section  we  discuss  a  generic  class  of  continuum  models  for 
flexible  structures  describing  the  abstract  mathematical  formulation  of  the 
models  as  a  framework  for  the  design  of  control  laws. 

In  the  second  section  we  show  how  direct  frequency  domain  designs  for 
control  la  vs  may  be  achieved  for  this  class  of  models  based  on  a  spectral 
factorization  procedure  which  replaces  the  usual  computation  of  Riccati 
equations.  This  procedure  uses  the  infinite  dimensional  (continuum)  model 
of  the  structure,  and  it  leads  to  distributed  (or  localized)  controls  which 
impact  the  macroscopic  behavior  of  the  structure  as  a  whole. 

In  the  third  section  we  examine  the  problem  of  deriving  transfer  function 
representations  of  the  structural  models  as  required  in  the  frequency  domain 
design  procedure.  We  show  that  the  analysis  of  certain  types  structures  - 
beams  with  one  or  more  degrees  of  freedom  -  can  be  “automated”  (though 
this  is  far  from  trivial)  using  symbolic  manipulation  systems  (SMP).  Hy¬ 
brid  systems  are  also  considered  and  models  are  developed  based  on  the 
interconnection  of  components  and  subsystems  with  a  careful  delineation 
of  the  causal  relationships  among  the  components. 

In  section  4  we  describe  an  analytical  procedure  for  the  derivation  of 
continuum  models  for  large  scale  structures  with  a  regular  infrastructure. 
We  focus  on  the  representation  of  dynamics  of  truss  systems  and  the  ther¬ 
mal  transport  properties  of  lattice  structures.  Control  problems  for  typical 
systems  of  this  type  are  also  discussed. 
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1  Introduction 


It  is  now  generally  accepted  that  large,  low  mass  density  lattice  structures 
will  be  essential  for  several  near  term  space  applications.  Moreover,  it  is 
apparent  that  active  control  of  structural  vibrations  will  be  necessary  to 
enhance  their  stiffness  and  damping  properties.  In  this  report  we  consider 
the  construction  of  effective  mathematical  models  for  elastic  dynamics  of 
space  structures  with  the  objective  of  designing  active  control  laws  for  these 
systems. 

The  success  of  active  control  for  such  structures  will  hinge  to  some 
extent  on  the  ability  of  a  control  law  to  react  to  vibratory  responses  which 
may  be  initially  localized  before  they  propagate  throughout  a  structure. 
This  leads  naturally  to  questions  of  how  to  implement  active  control  so 
as  to  distribute  the  control  effort  spatially  as  it  is  needed.  We  argue  that 
well  known  methods  exist  for  the  control  of  distributed  parameter  systems 
and  can  be  effectively  applied  if  continuum  models  for  the  candidate  space 
structures  can  be  computed.  The  nature  of  the  required  models  is  however 
quite  different  from  the  more  standard  finite  element  models  which  are 
popular  for  large  structural  analysis  problems  throughout  the  aerospace 
industry. 

We  begin  in  this  section  with  a  review  of  continuum  models  for  active 
structural  control.  We  highlight  the  nature  of  abstract  state  space  models 
for  these  systems.  In  this  study  we  have  employed  one  method  for  the 
computation  of  a  distributed  control  law  for  continuum  dynamics  based  on 
a  numerical  procedure  for  spectral  factorization.  This  method  requires  the 
computation  of  certain  representation  for  an  underlying  state-space  model 
for  the  structure  to  be  controlled.  In  a  later  section  we  discuss  the  theo¬ 
retical  basis  for  computing  effective  distributed  parameter  models  for  large 
truss  structures  with  regular  lattice  infastructure.  The  method  which  in¬ 
volves  “homogenization’’  (an  asymptotic  analysis  of  multiple  scales)  leads 
to  the  well  known  Timeshenko  model  for  beam  dynamics.  The  analysis  pro¬ 
vides  formulae  for  the  effective  beam  parameters  which  are  quite  different 
than  have  been  suggested  by  other  averaging  schemes  [20,21,22). 
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Comprehensive  models  of  flexible  spacecraft  dynamics  will  involve  sys¬ 
tems  with  fairly  complex  interconnections  of  lumped  and  distributed  sub¬ 
systems,  and  therefore,  we  intend  to  construct  the  overall  models  by  first 
developing  subsystem  models  and  then  combining  them  according  to  the 
required  interconnection  rules.  These  interconnections  lead  to  basic  ques¬ 
tions  of  causality  and  well-posedness  of  certain  standard  models  for  beams. 
These  questions  are  crucial  to  the  computation  of  hybrid,  state-space  mod¬ 
eling  of  an  integrated  space  platform. 

Throughout  this  effort  we  have  focused  on  the  potential  for  automatic 
and  computer-aided  computation  of  the  models  by  a  combination  of  modern 
computer  algebra  [23]  (symbolic  manipulation)  and  numerical  methods.  In 
our  efforts  we  have  used  the  program  SMP  [24,25].  We  will  review  the 
progress  in  using  SMP  for  the  computation  of  irrational  transfer  function 
models  for  hybrid  problems  in  a  later  section. 

1.1  Generic  Models  for  Structural  Dynamics 

In  this  section,  we  discuss  a  generic  model  for  elastic  dynamics  of  structures 
from  the  point  of  view  of  continuum  modeling.  We  will  summarize  the  con¬ 
struction  of  a  state  space  model  and  introduce  a  typical  control  problem 
for  vibration  suppression.  We  highlight  the  modal  approximations  which 
are  popular  for  these  problems  and  proceed  to  demonstrate  an  effective 
alternate  technique  for  model  construction  and  control  computation  based 
on  the  semi-group  property  [26]  of  a  state  space  model.  Effectively,  mod¬ 
eling  and  control  law  computation  can  proceed  in  the  frequency  domain, 
based  on  transfer  function  methods,  permitting  the  direct  computation  of 
a  resolvent  operator.  We  focus  on  the  class  of  structural  contro  1  problems 
for  which  the  question  of  control  of  propagation  of  wave-like  disturbances 
is  important.  In  this  framework,  we  can  present  the  semi-group  theory  by 
concrete  computations  of  practical  interest  to  structural  and  control  system 
engineers. 

The  standard  linear  continuum  model  for  a  flexible  structure  [27]  is 


(1) 


given  by  a  system  of  partial  differential  equations  (PDE) 

.  .d2w(t,z)  _  dw(t,z)  .  ,  \  . 

m («) — dt2  ;  +  -Do — +  A0w{t,z)  =  F(f,z) 

where  w[t,z)  is  an  N-vector  of  displacements  of  a  structure  0  with  respect 
to  some  equilibrium  for  fl  in  a  bounded,  open  set  in  ZN .  The  (vector) 
z  £  fl  a  coordinate  in  fl.  We  assume  the  boundary  dCl  is  smooth.  The  mass 
density  m(z)  is  positive  definite  and  bounded  on  dCl.  The  damping  term 
Dodtv/dt  contains  both  (asymmetric)  gyroscopic  and  (symmetric)  struc¬ 
tural  damping  effects.  The  internal  restoring  force  Aow  is  generated  by  a 
time-invariant,  differential  operator  A0  for  the  structure.  For  most  common 
structural  models,  Aq  is  an  unbounded  differential  operator  with  domain 
D(Ao)  consisting  of  certain  smooth  functions  satisfying  appropriate  bound¬ 
ary  conditions  on  dCl.  Thus,  for  these  problems,  D(A0)  is  typically  dense 
in  the  Hilbert  space  Ho  =  L2(Cl).  Often  (but  not  always),  the  spectrum  of 
A0 ,  a(A0),  consists  of  discrete  eigenvalues  with  associated  eigenfunctions 
which  constitute  a  basis  for  LJ(fl). 

The  applied  force  distribution  F(t,z)  can  be  thought  of  as  consisting  of 
three  components 

F{t,z)  =  Fd{t,z)  +  Fe(f,z)  -I-  Fa{t,z)  (2) 

where  Fd  is  N-vector  of  exogenous  disturbances  (possibly  forces  and  torques), 
Fe  is  a  continuous,  distributed  controlled  force  field  (an  available  option  in 
only  some  special  applications),  and  Fa  represents  controlled  forces  due  to 
localized  actuation; 

Fa{t,z)  =  £M2)uj(0  =  B0u(t).  (3) 

)=i 

The  actuator  influence  functions  bj(z)  are  highly  localized  in  fl  and  can  be 
approximated  by  delta  functions.  Measurements  are  available  from  a  finite 
number  p  of  sensors 

1/(0  =  C0w  +  (4) 

where  j/(l)  is  a  p-vector.  The  operators  Bo  :  — *  Mo,  C0  :  Ho  —>  kp,  and 

C'0  :  Hq  —*  kp  are  bounded. 


The  standard  optimal  regulator  control  problem  for  this  model  is  to  find 
the  controls  Uj(t),  j  =  (we  ignore  the  possibility  of  Fe)  given  the 

observations  y(t)  to  maintain  the  system  state,  e.g., 


such  that  the  performance  index 

J{u)  =  Jo  (N!S‘  +  euTu)dt  (6) 

is  minimized  where  ||i||q°  =  ( x,Qx)m  defined  in  an  appropriate  Hilbert 
space  for  the  state  x.  This  is  the  generic  control  problem  surveyed  in 
Balas  [28].  In  this  report,  we  will  concentrate  on  the  construction  of  Btate 
space  models  and  computational  aspects  of  equations  of  the  form  (l)  and 
of  optimal  (discrete)  controls  Uj(t)  appearing  in  (3). 


1.2  State  Space  Models 


The  choice  of  state  space  given  by  (5)  is  often  attractive  for  models  i  n 
the  generic  form  (l).  (We  will  discuss  later  that  attractive  alternate  state 
space  models  can  arise  in  hybrid  constructions.)  A  natural  assumption  for 
structural  problems  [27]  is  that  Ao  is  symmetric  with  compact  resolvent 
and  discrete  (real)  spectrum  which  is  bounded  from  below.  The  state  (5) 
can  be  considered  as  an  element  of  a  Hilbert  space  )(  =  D[A0 )x/2  x  with 
the  energy  norm 

IMl!  =  (xi,  A0xi)0  +  (mij.ij)  o  (7) 

where  the  first  term  represents  potential  energy  and  the  second  term  is 
kinetic  energy.  Thus  the  (abstract)  state  space  model  can  be  written 

dx(t,  z ) 


dt 


=  Ax{t,  z)  +  Bu(t) 


(8) 


where  y[t)  =  Cx(t,z) 


A  = 


0 

I 

,  B  = 

0  ' 

—Ao 

-Do, 

Bo, 

,  C  —  [Cq,  Cq 


(9) 


r 

lAjl 
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For  the  elastic  dynamics  of  space  structures,  there  is  always  some  (possibly 
small)  damping  D0  appearing  in  (1)  which  causes  A  to  b  e  dissipative.  Thus 
the  criteria  of  the  Hille-Yoshida-Phillips  theorem  [26]  are  satisfied  and  A 
generates  a  Co-semigroup  with  an  operator  which  we  write  suggestively  as 
eAt.  Moreover,  such  models  are  “hyperbolic”  [27]  in  the  sense  that  the 
semigroup  is  a  contraction,  i.e.,  ||ex‘j|  <  1  and  all  but  the  zero  frequency 
poles  are  only  slightly  damped;  ||eA‘||  <  t~st  for  some  small  6  >  0.  As 
we  shall  discuss  in  section  3,  some  popular  models  for  structural  elements 
such  as  beams  do  not  fit  in  this  framework.  However,  this  framework 
includes  models  appropriate  for  considerations  of  wave-like  dynamics  which 
propagate  causally  in  the  spatial  domain.  For  such  models,  the  question 
of  how  to  compute  controls  u(t)  and  system  response  x(f)  focuses  on  the 
so-called  mild  solution  of  (8); 

x(t,z)  =  eAtx(0,  z)  +  f  Bu{o)do.  (10) 

Jo 

Various  methods  are  available  for  approximation  of  the  system  (8)  Burns 
and  Cliff  [29].  One  popular  method  is  based  on  a  modal  (eigen-)  expansion 
of  A  which  generates  a  sequence  of  finite  dimensional  subspaces  >4  C  V. ,  k  — 
1,2,...,  where  Mt  =  span{^y,  j  =  1, . . . ,  fc}  and  the  <£;(z)  are  eigenfunctions 
(or  mode  shapes)  for  A.  Based  on  this  approximation,  a  sequence  of  finite 
dimensional  models  for  (5)  can  be  generated; 

i(*)(f)  =  A(*>x<*>(f)  +  B<*>u(t).  (11) 

Using  a  truncated  model  (11)  with  k  finite  and  the  performance  index 
J(u)  (6)  projected  onto  the  space  #*,  one  can  solve  the  associated  optimal 
control  problem  for  the  first  k  modes  of  A.  However,  as  noted  in  Balas 
[28]  in  all  but  a  few  special  cases,  the  control  law  when  applied  to  the 
system  (5)  will  excite  higher  order  modes.  The  inherent  robustness  and 
stability  properties  as  well  as  the  degree  of  suboptimality  of  control  laws 
based  on  such  truncated  modal  approximations  has  received  a  great  deal 
of  attention  in  both  the  engineering  and  mathematics  literature  [27,28]. 
Various  alternate  approaches  are  available  which  deal  directly  with  infinite 
dimensional  control  problem  given  by  (6)  and  (10)  -  at  least  abstractly 


"I 


I 

N 

*> 

I 


P 


Cn 


Russell  [30].  One  method  suggested  by  Davis  [31,32]  offers  the  advantage 
of  a  computational  procedure  for  approximating  the  true  optimal  control 
in  terms  of  the  required  control  bandwidth.  The  method  is  based  on  an 
extension  of  a  Wiener-Hopf  solution  [32]  for  the  abstract  control  problem. 


2  Wiener-Hopf  Methods  for  Computation  of 
Optimal,  State-Feedback  Control:  The  Davis 
Stenger  Algorithm 


The  connections  between  least  squares  optimization,  spectral  factorization, 
and  algebraic  Riccati  equations  have  been  considered  important  in  control 
theory  for  many  years.  (See  e.g.,  Anderson  [66],  Brockett  [33],  Willems  [69], 
Helton  [68],  and  the  references  therein).  To  see  how  the  connection  arises, 
consider  the  standard,  finite-dimensional,  infinite  time  regulator  problem: 

min  /“[Moil2  +  ||y(OI|a«&  (1) 

uGli**  JO 

subject  to  the  linear,  time-invariant  system  model 

x(t)  =  Ax(t)  +  Bu(t),  x(0)  =  x0  (2) 

and  controlled  output 

y(t)  =  Cx(t),  t>  0.  (3) 

The  transfer  function  relating  the  Laplace  transform  of  the  input  vector 
u(s)  to  the  output  vector  y(s)  is 

y[s)  =  G(s)u{s), 

G{s )  =  C\sI-A]~lB.  (4) 

The  optimal  control  is  known  to  be  a  linear  state  feedback 

u(t)  =  -Kopix(t)  =  -BTPx{t)  (5) 

where  P  is  the  unique,  positive  definite  symmetric  solution  to  algebraic 
(matrix)  Riccati  equation, 

PA  +  AtP  -  PBtBP  +  CTC  =  0.  (6) 

Standard  algebraic  manipulations  based  on  (l)-(6)  provide  the  spectral 
factorization  relation 


[/  +  K^i-al-  A)-lB}T\I+Kopt(aI-  A)~l B\ 
=  I  +  Bt(-s  -  At)~*CtC(sI  -  A)~lB. 


[57,  pp.  68].  Clearly  (7)  cam  be  rewritten 

H(s)  =  I  +  Gt{-s)G{s)  =  Ft{-s)F{s )  (8) 

where  .F(s)  =  I  +  Kopt[sI  —  A]~1B  is  the  causal  spectral  factor  of  H(s). 
In  fact  many  desirable  properties  of  linear  quadratic  regulator  design  (e.g. 
inherent  robustness  properties)  follow  from  (8)  and  the  interpretation  of 
F(s)  as  the  optimal  return-difference  operator  for  the  problem  (l)-(3). 

Recall  that  for  a  closed  loop  control  given  by  the  transfer  function  re¬ 
lations 

y{s)  =  Q(s)u(s) 

u{s)  =  f(s)  -  e(s)  (9) 

that  the  return-difference  with  respect  to  the  returned  signal  y(s)  —  e(s) 
is  computed  as  the  difference  between  y(s)  and  u(s)  if  the  loop  is  broken 
there.  Thus, 

SW  -  »(«)  =  [X  +  CMI  f(»). 

In  (8)  F(s)  is  the  optimal  return-difference  for  loop  breaking  at  the 

control  with  state  feedback.  We  remark  that  for  continuum  models  for 
flexible  structures  that  the  underlying  state  space  models  are  infinite  dimen¬ 
sional  and  F(s)  represents  an  ideal  limit  on  achievable  control  performance 
with  state  feedback. 

In  Nyquist  stability  theory  the  complex  contour  1  +  H(ju)  for  u  G  S? 
is  called  the  Nyquist  contour.  Practical  aspects  of  Nyquist  stability  tests 
have  been  the  cornerstone  of  control  system  design  for  at  least  50  years. 
These  stability  tests  have  been  extended  to  include  an  important  class  of 
distributed  parameter  systems  where  the  resulting  return  difference  is  an 
irrational  transfer  function  [58].  For  our  purposes  in  this  study  the  Nyquist 
test  will  prove  central  and  we  will  focus  on  the  class  of  irrational  transfer 
functions  for  which:  (l)  spectral  factorization  can  be  effectively  computed 
at  a  finite  number  of  samples  and  (2)  the  Nyquist  theory  of  stability  is  well 
defined. 


In  addition  to  (7)  a  useful  integral  equation  can  be  derived  for  the 
optimal  state  feedback  under  the  additional  assumption  that  the  spectrum 
of  the  operator  A  in  (2)  is  contained  the  open  left  half  plane  (C_).  We  will 
state  the  result  as  a  theorem  and  outline  a  proof  for  the  finite  dimensional 
case.  Next  we  consider  the  additional  assumptions  necessary  for  the  result 
to  hold  for  control  of  a  distributed  parameter  problem  in  the  form  of  (1)- 

(3)- 

Theorem  1  (Davis)  With  the  optimal  linear  regulator  problem  as  defined 
in  (l)-(S)  the  optimal  feedback  control  (if  it  exists)  u(t)  =  Koptx(t),  can  be 
computed  as 

Kopt  =  ~  l°°  [F*(iw)]_1  G*{iu>)CR{iu>,  A)du.  (10) 
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under  the  assumption  that  A  is  a  stable  (matrix)  operator. 

In  (10)  we  use  the  notation  G(s)  as  defined  in  (4),  F*(iu)  —  FT(— tu;)  with 
.F(s)  the  causal  spectral  factor  given  by  (8)  and  P(iu>;  A)  =  [twJ  —  A]~l  the 
(matrix)  resolvent  for  A. 

Proof:  From  standard  results  [33]  the  optimal  control  (if  it  exists)  will 
stabilize  the  closed  loop  system  so  that  a(A  —  BKopt)  C  C_  where  C_ 
is  the  open  left  half  of  the  complex  plane.  Construct  a  closed  rectifiable 
contour  T  in  the  complex  plane  consisting  of  a  relatively  large  portion  of 
the  *w-axis  and  a  semicircular  portion  in  the  left  half  plane  such  that  T 
encircles  a(A  —  BKopt )  in  the  positive  sense.  Let  Ae  =  A  —  BK^f  Then 

-  A,r  = i.  (in 

By  assumption  o(—AT)  is  contained  in  C+  and 

^^[s/+At]  'ds  =  0.  (12) 


From  (6)  we  write 


ATP  +  PAe  =  -CTC. 


(13) 


>■1 


This  leads  to  the  relation 


P{sl  -  Ac)"1  +  {-si  -  AT)~lP  =  {-si  -  At)~1CtC{sI  -  A,)'1-  (14) 
Now  integrate  (14)  on  T  and  use  (11)  and  (12)  to  get 

p  =  jlj  jf  (-*/  -  AT)-'CTC(.I  -  A,)-'ds.  (15) 

Since  the  optimal  gain  is  Kopt  =  BTP  =  \PB\T  we  get 

Kopt  =  ^-}T  BT{sI  ~  Ae)~lCTC{—sI  -  AT)~1ds.  (16) 


Now  from  (7)  and  (8)  we  get  that 

C{sl  -  At)~lB  =  C{sl  -  A)_1J5  [/  +  Kopt{sI  -  A)_1JB] 


and  therefore 


Thus  we  can  write 


C{sI-Ae)~1B  =  G(s)F-x(s). 


Kopt  =  2 hi  F~T(a)GT(*)C(-aI  ~  ^)_1^.  (18) 

Finally  (10)  is  determined  by  substituting  s  =  iu;  under  the  observation 
that  G{iu)  -tOuw-*  oo. 

For  (10)  to  hold  when  G{s)  is  irrational  we  must  impose  some  constraints 
on  the  spectral  properties  of  the  associated  infinite-dimensional  operator  A. 
Consideration  for  the  computation  of  the  Riccati  operator  P  by  the  integral 
formula  (15)  suggests  that  the  spectrum  of  A  must  consist  of  a  countably 
infinite  set  of  eigenvalues  so  that  the  path  integrals  can  be  computed  as 
in  (11)  and  (12).  In  addition,  the  spectral  sets  o{At )  and  c{A )  must  be 
separated  by  the  imaginary  axis  (including  the  point  at  infinity).  Also, 
the  observation  that  G(»w)  — ►  0  as  u  — ►  oo  for  rational  functions  must  be 
replaced  by  an  equivalent  assumption  which  further  restricts  the  class  of 
irrational  transfer  functions. 


Wi 
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Effectively  these  additional  assumptions  will  be  guaranteed  by  the  class 
of  transfer  functions  for  which  spectral  factorization  can  be  performed. 
We  will  consider  spectral  factorization  for  distributed  systems  next.  More 
importantly  with  the  usual  assumptions  used  in  constructing  continuum 
models  for  mechanical  structures  it  appears  that  the  resulting  transfer  func¬ 
tions  will  have  the  appropriate  properties.  However,  in  this  study  we  have 
encountered  much  confusion  in  the  literature.  As  discussed  in  the  intro¬ 
duction  a  major  goal  of  this  study  was  to  provide  a  consistent  method  for 
model  construction  leading  to  an  appropriate  class  of  transfer  functions  for 
models  which  are  both  well-posed  and  have  appropriate  spectral  proper¬ 
ties.  In  the  next  section  we  delineate  the  specific  class  of  transfer  functions 
for  which  spectral  factorization  can  be  computed  efficiently  by  a  numerical 
algorithm.  We  indicate  the  basis  for  the  algorithm  and  discuss  the  im¬ 
plications  of  sampling  and  interpolation  of  the  spectral  factor.  Finally  we 
discuss  the  relationship  between  rational  approximation  and  modal  control. 


2.1  Spectral  Factorization 

In  this  section  we  review  the  basis  for  an  interactive  algorithm  for  compu¬ 
tation  of  a  frequency  sampled  spectral  factor.  The  algorithm  (due  to  Davis 
and  Dickinson  [32])  provides  an  effective  computational  tool  for  obtaining 
the  optimal  gain  Kopt  via  (10)  without  regard  to  computational  difficulties 
associated  with  large  dimensional  Riccati  equations.  Convergence  of  the 
algorithm  depends  on  certain  technical  assumptions  which  delineate  the 
class  of  transfer  functions.  In  the  finite  dimensional  setting,  a  recursive 
algorithm  for  computation  of  the  causal  spectral  factor  F(s)  follows  from 
a  Newton-Raphson  iteration  for  the  matrix  Riccati  equation  (6)  given  an 
initial  stabilizing  feedback  K0  =  B*P0\ 

Pn+l(A  -  BB*P„)  +  (A  —  BB*PnyPn+1  =  —C*C  -  PnBB*Pn. 


At  the  nth  iteration  one  can  take  an  approximation  to  the  causal  spectral 
factor  as 


Then  following  Davis  and  Dickinson  [32]  this  leads  to  the  form  of  the  algo¬ 
rithm 

f.+>M  =  K  {|JCM|-,H(.'u.)|F„(.u,)]-1}  *•„(«,),  (19) 

where  P+  is  the  causal  projection  operator  defined  on  the  convolution  alge¬ 
bra  I  ©  L\  or  I  ©  L2  by 

P+  {/  +  f{t)e~iutdt}  =  /  +  jH  f(t)e~iutdt. 

Computation  of  causal  projection  of  a  signal  is  a  standard  problem  in 
signal  processing  which  can  be  effectively  solved  through  the  use  of  a  Hilbert 
transform  [59],  [60]. 


Definition  .1  The  Hilbert  transform  of  a  signal  f(t)  is  given  as 


/(*0 
t  —  T 


dr  =  /(f) 


* 


l_ 
rt  ’ 


(20) 


The  basic  utility  of  the  Hilbert  transform  is  evident  by  examination  of  its 
Fourier  transform. 


Fact  2  Let  u)  be  the  Fourier  transform  of  }(t)  and  F(u>)  be  the  Fourier 
transform  of  f(t). 

p(u)  =  r  -  r  ji^-dr  eiutdt.  (21) 

J- oo  T  J- 00  t  —  T 


Since  for  H(u)  =  ±etwtdt,  we  see  that  H(u)  is  a  complex  function  with 
the  properties 


(23) 


I 


4* 


Therefore  from  (21)  and  (22)  we  get  that 


bluj)  _  /  -«^(w)  w  >  0 
“  {  iF(u)  u  <  0 


Now  to  compute  the  causal  projection  of  f(t)  (as  in  (19))  given  frequency 
domain  data  F(u>)  we  first  compute  the  Hilbert  transform 

P(u)  =  F(u>)  *  ~ 

7UO 

By  duality  of  the  Fourier  transform  pair  the  property  (23)  holds  in  th  e 
domain  t\ 

7~l  {^(w)}  =  — isgn (0/(0* 

Finally  causal  projection  can  be  computed  as 


^{F(«)>-l[F(w) +  <!•(«)]. 


Before  we  consider  computational  issues  further  we  review  the  extension  of 
this  algorithm  to  irrational  transfer  functions.  The  questions  of  existence 
and  uniqueness  of  the  spectral  factorization  of  the  transform  H(s)  =  I  -+- 
GT(s)G(s)  naturally  lead  to  conditions  for  which  S(u>)  =  GT(— iw)G(iu)  is 
positive  semidefinite  for  lj  real.  In  the  convergence  proof  of  the  iteration 
(19)  Davis  assumes  that  G(s)  is  the  transform  of  a  rea  1,  vector- valued 
function  which  is  both  integrable  and  square  integrable;  i.e.,  G( »w)  G  7(LiCi 

Lj). 

With  these  assumptions  it  is  clear  from  the  classical  theory  of  Gohberg 
and  Krein  [54]  that  H(s)  has  a  unique  spectral  factorization  as  given  in  (8) 
with 

F±{iu)  -  I  e  7{Lt) 

where  7{L>i)  is  the  class  of  Fourier  transforms  of  functions  in  L\  with 
positive  support.  As  noted  in  [32],  the  assumptions  on  G(<)  in  fact  imply 
that 

F±1(tu;)  —  /  G  J(L+nI+) 

and  F(iui)  =  F(— iw). 
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Using  this  assumption  Davis  is  able  to  show  that  the  recursion  (19) 
starting  from  an  initial  F0(iu)  €  /(LinZ^)  has  all  iterates  Fn(»w)  €  7{L\  n 
L2)  and  that: 


limn— » oe  F*(w)F„(u;)  =  H{  u)  almost  everywhere 
limn^oo  Fn(s)  =  F(s)  for  all  s  with  SRe  s  >  0. 


(25) 

(26) 


These  results  should  not  be  surprising  for  the  class  of  transfer  functions 
considered.  The  following  theorems  summarize  well  known  properties  of 
Li  and  L2  functions. 


Theorem  3  ([68])  If  f  e  Lx  then 


1.  cj  *— ►  /(u>)  is  uniformly  continuous  for  u>  E  R 


2.  /(<*>)  €  L0 


S.  |/(w)|  —*  0  as  |w|  — ►  00 

4.  f{t)  —  A  f(iu)eiwtdu  almost  everywhere  in  R. 


1 


i 
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So  we  see  that  /  €  L\  means  that  f  (*u>)  is  bounded  on  u;  and  therefore  the 
Fourier  transform  is  well  defined  while  /  £  Lj  provides  consistent  approxi¬ 
mation  theory  for  band  limited  signals.  We  note  that  the  third  Li  property 
means  that  such  transfer  functions  are  effectively  band  limited  and  strictly 
proper. 

Application  of  spectral  factor  to  the  integral  formula  (10)  will  further 
restrict  consideration  to  G(s)  both  causal  and  stable  so  that  G(s)  is  analytic 
for  3?es  >  0.  Thus  the  transfer  functions  we  are  considering  belong  to  the 
Hardy  space  G(s)  €  H2  n  H00.  Recall  that  by  definition  f  6  H2  if  /  is 
complex  valued  and  analytic  in  C+  (the  open  right  half  of  the  complex 
plane)  and 

/oo  „ 

\n°  +  iu)\du>  <  oo. 

-oo 

while  /  €  H°°  if  f  is  complex  valued  and  analytic  in  C+  and 

sup  \f{s)\  <  oo. 

«eC+ 

The  first  property  is  inherited  from  /  €  Lj  while  the  second  comes  from 
/  €  Li. 


Finally,  we  remark  that  for  certain  mechanical  structures  that  the  trans¬ 
fer  function  models  can  be  factorized  in  a  “product  expansion”  [7l] 


/w = n 


*=i 


1  +  (fVfg) 

i  +  («2M)' 


(27) 


Thus  we  conclude  these  remarks  by  indicating  that  the  class  of  transfer 
function  models  for  which  spectral  factorization  can  be  computed  by  the 
present  method  are  meromorphic  and  have  representations  as  (27). 

2.1.1  Remarks  on  Algorithm  Construction 


With  regard  to  the  integral  formula  for  the  optimal  state  feedback  gain  it 
is  clear  that  we  need  [F(iu;)]_1.  Thus  as  suggested  by  Davis  it  is  convenient 


to  implement  the  iteration  in  the  form 


|f»+.r'  =  [FJ-1  (7  +  F+  {[F„-l-  H [F„|-‘  -  /})’■ .  (28) 

Furthermore  by  initializing  with  Fq  a  diagonal  matrix  with  diagonal  ele¬ 
ments  equal  to  the  spectral  factors  of  the  diagonal  elements  of  H  the  Becond 
term  of  (28)  remains  a  perturbation  of  the  identity  (since  [F*]-1.ff[Fn]-1  — 
/  — ♦  0)  which  regularizes  the  computations.  The  diagonal  initialization 
guarantees  that  the  first  residual  [Fo]-1fl’[Fo]-1  has  ones  on  the  diagonal 
and  all  off  diagonal  elements  less  than  one  in  magnitude.  Using  the  proper¬ 
ties  of  the  Hilbert  transform  and  the  formula  (24)  one  can  readily  compute 
the  causal  spectral  factor  for  the  individual  scaler  transfer  functions  di¬ 
rectly  (i.e.  without  iteration).  In  particular,  the  kth  diagonal  element  of 
H(u),hk(u)  is  a  real  valued  function  with  hk(ui)  >  0.  Let  hk( w)  be  the 
Hilbert  transform  of  hk(u)]  viz., 

hk{u)  =  -  r  da  (29) 

ir  J-oo  u>  -  a  v  ' 

then  the  causal  spectral  factor  hk( iw)  =  /*(— iu)fk{iu})  is  given  by 

fk  M  =  y/Xe  Mw)  e"  ^  A(w)  •  (30) 


Finally,  we  remark  that  numerical  computation  of  the  Hilbert  trans¬ 
form  can  be  achieved  in  several  ways.  Direct  numerical  integration  of  (29) 
is  complicated  by  the  fact  that  the  integral  is  convergent  in  the  Cauchy 
principal  value  sense.  Effective  quadrature  algorithms  for  such  problems 
have  been  coded  and  tested.  A  public  domain  version  utilizing  an  adap¬ 
tive  quadrature  procedure  is  available  in  the  routine  QAWC  contained  in 
a  software  package  called  QUADPACK  J61). 

Another  approach  is  taken  by  Davis  [32]  based  on  an  algorithm  of 
Stenger  [62].  This  procedure  essentially  implements  a  discrete  (sampled) 
version  of  the  required  computation  using  a  digital  Hilbert  transform.  For 
a  finite  number  of  sample  points  the  Hilbert  transform  computation  can 
be  implemented  by  taking  a  discrete  “fast  fourier  transform”  (FFT)  of 
the  sampled  data  and  shifting  the  imaginary  part  of  the  transformed  data 


according  to  (21).  It  is  well  known  that  control  of  error  induced  by  the 
sampling  process  (Gibbs  phenomenon)  requires  the  careful  choice  of  “data 
windows"  or  weighting  functions  for  the  computation.  Although  these  con¬ 
siderations  are  well  described  in  the  literature  on  digital  signal  processing 
[60],  [64],  [63]  it  is  not  apparent  that  these  considerations  were  contained 
in  the  work  of  Stenger. 

In  our  experience,  direct  implementation  of  the  algorithm  described  in 
Davis  and  Dickinson,  based  on  the  causal  projection  of  F.  Stenger  [62]  tends 
to  be  unreliable.  Since  detailed  design  of  window  functions  for  discrete 
Hilbert  transforms  was  considered  outside  the  scope  of  the  present  study 
we  have  found  it  expedient  to  use  the  adaptive  quadrature  software  from 
QUADPACK  [61]. 

2.1.2  Relationship  to  Modal  Methods 

In  applications  one  typically  sees  the  use  of  modal  methods  in  the  modeling 
and  control  of  continuum  models  for  structures  with  a  procedure  as  follows 
[28], [72]: 

1.  Determine  a  finite-dimensional,  modal  (or  finite  element)  approxima¬ 
tion  to  the  distributed  parameter  problem  yielding 

x(<)  =  Ax(t)  +  Bu(t) 

with  z(t)  a  finite-dimensional,  state  vector  representing  displacements 
and  velocities  for  a  finite  number  of  structural  modes. 

2.  The  finite-dimensional  model  may  be  further  truncated  to  achieve  a 
“reduced-order”  approximation  which  will  be  used  for  control  synthe¬ 
sis. 

3.  Standard  linear  quadratic  regulation  theory  is  applied  to  the  reduced- 
order  model.  The  resulting  matrix  Riccati  equation  can  be  solved 
(in  principal)  although  at  some  significant  computational  cost  for  an 
arbitrary  (large)  number  of  modes. 


4.  Implementation  of  the  active  control  law  requires  on-line  estimation  of 
the  particular  plant  modes  included  for  control  design  A  subsequent 
analysis  of  the  effect  of  control  action  on  the  residual  (truncated) 
model  may  force  a  rather  drastic  reduction  in  the  effective  loop  gain 
to  avoid  a  potentially  destabilizing  phenomenon  known  as  “spillover” . 


We  remark  that  control  laws  of  this  type  involve  state  feedback  for 
a  truncated  model  and  therefore  require  real-time  measurement  (or  esti¬ 
mation)  of  certain  individual  modes.  Since  both  mode  shapes  and  modal 
frequencies  can  change  under  the  influence  of  feedback  control  such  feed¬ 
back  can  be  difficult  to  realize.  This  is  further  complicated  by  the  fact  that 
typically,  exact  modal  data  is  not  available  in  the  construction  of  step  1 
using  finite  element  procedures.  The  Wiener-Hopf  methods  we  are  inves¬ 
tigating  here  are  based  on  quite  different  modeling  schemes  and  we  will 
discuss  modeling  in  more  detail  in  later  sections.  However,  since  it  is  often 
desired  to  control  “certain  troublesome  modes”  in  engineering  design  it  is 
useful  to  consider  the  relationship  of  these  methods  with  the  objective  of 
modal  control. 

Consider  the  state  space  continuum  model  for  a  structure  considered  in 
section  1. 

— =  A  x[t,  z)  +  Bu(t)  (31) 

where  x(t,z)  6  )/(0),  a  real  Hilbert  space  satisfying  boundary  conditions 
for  z  6  fi  is  finite-dimensional  control,  A,B  are  operators  on  W(fi)  and  we 
assume  that  for  structures  A  is  densely  defined  on  M  and  maximally  dissi¬ 
pative  with  a  discrete  spectrum  consisting  of  its  eigenvalues.  The  solution 
to  (31)  can  therefore  be  given  as 

x(t,  z)  =  eAfz(0,  z)  +  f  u(s)ds  (32) 

Jo 

where  eAt  is  the  semigroup  operator.  Roughly,  since  A  is  assumed  to  have 
discrete  spectrum  a  modal  analysis  can  be  based  on  an  eigen  expansion  of 
A  in  terms  of  eigenvalues  A*  and  eigenfunctions  <Pk(z) 


A  <pt{z)  =  A  kVkiz)- 


Considering  the  Laplace  transform  of  (32) 

x(s,z)  =  R(s\  A)x(0,  z)  +  ffBC(s,z)u(s)  (33) 

=  J  Gr(s,  z,  tu)x(0,  w)dw  +  Hbc{s,z)u(s),  (34) 

where  J?(s;A)  is  the  resolvent  for  A  which  can  be  given  by  the  integral 
operator  (as  shown)  with  kernel  called  the  Green’s  function. 

In  the  fram'"vork  of  the  distributed  parameter  (DP)  control  problem 
“modal  control”  as  an  objective  means  that  we  choose  the  output  for  control 
(2)  by  defining  the  projection  nN  :  M  HN  c  U  where  is  a  finite 
dimensional  “modal”  subspace  of  H.  Then  taking 

y(t)  =  CnNx(t,z) 


leads  to 

y(s,z )  =  j  irffGr(s,z,  w)x(0,  w)dw  -f  nNHBc{s,z)u{s)f 

where  the  elements  of  the  matrix  Green’s  function  can  be  given  in  terms  of 
a  corresponding  modal  expansion; 


nr-  ( . *M 

I't’r  Z,W )  —  2^  7~T 

M«) 


*=1 


with  p*(s)  =  sJ  +  2f*w*s  +  w*.  Here  are  the  frequency  and  damping 

of  the  eigenmodes  and  <Pk[z)  are  the  eigenfunctions  (modeshapes). 


To  compute  the  optimal  control  by  Wiener-Hopf  methods  we  must  per¬ 
form  spectral  factorization  on 


H{s)  =  I  +  Gt(-s)G{s) 

where  G(s)  =  CHbc[s,  z)  so  that  G(a)  is  a  rational  transfer  function. 
Computation  of  the  optimal  state  feedback  via  (10  requires  the  resolvent 

CR(iu>,  A)  z(0,  w)  =  J^7rNGr{8,z,vu)x(0,w)dw 
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which  is  also  rational  in  the  Laplace  variable  s.  Under  the  assumption  that 
all  truncated  modes  are  stable  then  the  resulting  control  will  stabilize  the 
system.  In  this  case  the  Davis  method  produces  an  effective  distributed 
feedback  control 

u(t)  =  Koptx[t,z)  =  [  K(z)x(t,z)dz 

for  the  truncated  modal  subspace  which  is  effectively  the  same  control  one 
would  get  by  the  procedure  described  earlier  assuming  that  exact  modal 
data  is  available.  However  notice  that  the  control  is  a  full  (distributed)  state 
feedback  so  that  (in  principle)  modal  measurements  are  not  required.  Fur¬ 
thermore,  computation  of  large  dimensional  Riccati  equations  is  replaced 
by  the  procedure  for  spectral  factorization.  We  remark  that  a  variety  of 
algorithms  exist  for  spectral  factorization  of  rational  functions  [70],  [73]. 

The  point  to  focus  on  is  that  the  method  of  Davis  completely  separates 
the  question  of  rational  approximation  (which  may  be  based  on  modal 
expansion)  of  the  model  from  the  computational  issues  associated  with 
determining  the  control.  Instead,  approximation  is  based  on  sampling  and 
interpolation  of  the  transfer  functions  is  involved. 

2.2  Importance  of  Damping  in  Models  for  Distributed 
Control 

For  the  present  study  we  have  attempted  to  compute  distributed  controls 
for  several  examples.  Several  negative  results  lead  us  to  reevaluate  the 
theoretical  basis  for  control  and  modeling  problem.  Our  difficulties  stem 
from  precise  computation  of  the  irrational  transfer  functions  for  several 
distributed  models  which  will  be  discussed  later.  Thus  our  control  laws  are 
computed  and  tested  on  bona  fide  distributed  parameter  models  and  not  on 
finite  dimensional  approximations;  therefore  the  idiosyncracies  of  various 
(academic)  examples  will  lead  to  obvious  difficulties.  In  this  section  we 
review  our  most  recent  conclusions  about  the  assumptions  in  modeling  for 
mechanical  structures  and  resulting  computation  of  optimal  control  laws. 


The  most  pertinent  comments  in  the  literature  appear  to  be  summa¬ 
rized  in  the  work  of  Gibson  [65].  In  particular,  the  importance  of  damping 
and  the  ability  to  compute  distributed  control  laws  for  infinite  dimensional 
models  for  mechanical  structures.  It  is  true  that  computation  of  control 
laws  for  such  problems  will  always  involve  approximation  in  modeling,  con¬ 
trol  problem  formulation,  and  numerical  computation.  Gibson  focuses  on 
this  question  of  approximation  in  modeling  and  control  for  distributed  pa¬ 
rameter  problems  via  the  use  of  finite  dimensional  models  and  control  com¬ 
putation  based  on  these  models.  Although  on  the  surface  the  Davis  method 
(based  on  spectral  factorization  of  irrational  transfer  functions)  is  not  sub¬ 
ject  to  finite  dimensional  approximation  (at  least  in  the  sense  of  modal 
expansions)  it  is  clear  that  approximation  via  sampling  and  interpolation 
of  the  transfer  functions  involved  is  required. 

Let  us  consider  some  observations  of  Gibson  [65].  First,  Consider  the 
class  of  elastic  mechanical  structures  modeled, 

z(t,  z )  +  C0i(t,  z)  +  A0x(t,z)  =  Bu(t)  (35) 

where  x  G  "H  (H)  appropriately  chosen  to  match  required  boundary  condi¬ 
tions,  u(t)  is  a  finite  dimensional  vector  of  controls,  A0  is  a  self  adjoint 
operator  Aq  :  D(A0)  — ♦  M  densely  defined  on  U.  Gibson  further  assumes 
that  Ao  is  coercive ;  i.e. 

(A0x,x)#  >  p2||x||#,  x  G  D(A0) 

and  Aq1  is  compact.  Co  is  nonnegative,  symmetric  linear  operator  and 
there  exists  -y  >  0  such  that 

||C0x||  <  7J||A0x||#,  x  G  D(A0).  (36) 

Finally,  Bo  is  taken  to  be  a  bounded  operator.  Gibson  shows  that  (36)  is  a 
necessary  condition  for  the  ~esulting  semigroup  operator  for 

A  =  [  0  1 

-A0  -Co  . 

on  M  x  H  to  be  uniformly  exponentially  stable4,  i.e.,  there  exists  M  >  0,  a  >  0 
such  that  ||e'4f||  <  Me~at  for  t  >  0.  Such  behavior  corresponds  to  the  case 
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when  damping  provides  a  uniform  decay  rate.  Gibson  addresses  the  ques¬ 
tion  of  convergence  and  stability  of  a  sequence  of  finite  dimensional  (pos¬ 
sibly  modal)  approximate  control  problems  to  the  unique  optimal  control 
for  the  distributed  model.  His  results  indicate  [65,  thm  4.1]  that  for  the 
quadratic  linear  regulator  problem  for  the  distributed  model  above  (35)  that 
if  a  solution  exists  the  exact  optimal  feedback  control  provides  a  closed  loop 
system  whose  infinitesimal  generator  Aci  =  A  +  BKopt  generates  a  strongly 
continuous  semigroup  which  is  uniformly  exponentially  stable. 

For  systems  without  damping;  i.e.  C0  =  0  in  (35),  Gibson  shows  that 
there  can  be  no  nonnegative,  selfadjoint  solution  for  the  algebraic  Riccati 
equation  for  the  distributed  model.  This  follows  from  the  observation  [65, 
theorem  5.13]  that  for  C0  =  0,  the  semigroup  generated  by  A  +  £  for  any 
£  a  compact  linear  operator,  cannot  be  uniformly  exponentially  stable. 

For  systems  with  damping  Gibson  shows  that  the  damping  must  be 
such  that  A  (the  generator  for  the  open  loop  system)  is  uniformly  exponen¬ 
tially  stable  in  order  that  uniform  exponential  stability  of  the  closed-loop 
system  can  be  obtained  by  compact  linear  feedback.  The  significance  of 
this  fact  is  that  the  physical  nature  of  damping  must  be  considered  in  the 
computation  of  distributed  parameter  control  and  for  problems  where  the 
available  control  affectors  provide  ‘localized’  (in  the  spatial  domain)  forces 
(or  torques).  Damping  with  a  uniform  decay  rate  (for  essentially  all  modes) 
is  required  for  such  systems  to  be  controlled  with  uniform  exponential  sta¬ 
bility.  Although  it  is  generally  agreed  that  physical  structures  will  have 
this  property  is  somewhat  disconcerting  that  many  standard  models  for 
simple  structural  elements  like  beams  with  internal  damping  do  not  have 
these  properties.  In  fact  with  the  exception  of  viscous  damping  there  does 
not  appear  to  be  a  single,  universally  excepted,  well-defined  mathematical 
model  for  beam  dynamics  which  is  obviously  appropriate  for  the  study  of 
distributed  parameter  control  of  structures. 

In  a  later  section  we  will  discuss  damping  models  for  structural  elements 
in  detail.  At  this  point  we  remark  that  in  this  study  it  became  apparent 
that  the  performance  of  active  control  of  elastic  structures  was  seen  to  be 
heavily  dependent  on  the  type  of  damping  models  used.  If  one  is  willing  to 
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assume  a  sufficient  amount  of  viscous  damping  then  stabilizing  control  can 
be  computed.  Since  viscous  effects  will  not  be  present  in  space  applications 
one  is  faced  with  a  choice  of  several  models  for  internal  damping.  Many 
of  these  models  lead  to  problems  for  which  either  a  uniform  decay  rate  is 
not  available  or  for  which  the  spectrum  of  the  operator  contains  more  than 
mere  eigenvalues.  These  models  cannot  be  stabilized  by  compact  linear 
feedback  whether  it  is  computed  by  Weiner-Hopf  methods  or  by  modal 
approximation  and  solution  of  a  finite-dimensional  quadratic  regulator  for 
the  reduced  model.  However  if  one  takes  the  second  path  the  relevant 
stability  questions  are  completely  lost  in  the  model  reduction. 

What  is  an  appropriate  choice  for  internal  damping  in  space  structures 
is  definitely  an  open  question.  This  issue  has  apparently  not  received  a 
great  deal  of  attention  in  the  aerospace  industry  especially  in  that  portion 
of  the  community  involved  in  control  of  large  flexible  structures.  These 
standard  procedure  here  (as  evidenced  in  for  instance  the  ACOSS  program) 
is  to  assume  “low”  modal  damping  can  be  added  to  the  reduced-order 
model  prior  to  control  system  design.  In  the  present  frequency  domain 
computations  one  is  forced  to  resolve  these  issues  before  proceeding  with 
control  computation. 


3  State  Space  Models  for  Distributed  Structural 
Elements 

3.1  Standard  Forms  for  Linear  PDEs 


In  this  section,  we  will  consider  the  problem  of  deriving  the  transfer  matrix 
description  for  typical  distributed  elements.  It  is  our  contention  that  the 
systems  of  interest  to  us-specifically  beams  with  one  space  variable  and 
perhaps  several  degrees  of  freedom-can  be  represented  by  one  of  two  stan¬ 
dard  forms.  Once  identifying  the  structure  of  these  standard  models,  it  is 
straightforward,  although  far  from  trivial,  to  mechanize  the  construction 
of  the  required  transfer  matrices  using  symbolic  computation.  Moreover, 
in  order  to  assemble  hybrid  system  models  by  the  interconnection  of  com¬ 
ponents  or  subsystems,  it  is  essential  to  have  a  clear  understanding  of  the 
causal  requirements  of  the  component  mathematical  models.  The  follow¬ 
ing  paragraphs  develop  the  required  concepts  in  terms  of  commonly  used 
structural  elements.  Since  typical  elements  interact  at  physical  boundaries, 
our  foremost  concern  is  with  the  formulation  of  appropriate  boundary  con¬ 
ditions  for  well-posed,  initial-boundary  value  problems. 


Before  proceeding,  we  establish  some  basic  terminology  associated  with 
systems  of  partial  differential  equations.  Consider  the  system  of  first  order 
partial  differential  equations  defined  for  t  >  0  and  0  <  x  <  L 


^dw  A .  dw  „„ 

E—  =  A  —  +  BV,  w  £  R 
at  dx 

If  E  is  nonsingular,  then  (l)  can  be  written 

dw  dw 

- —  —  A - (-  Bw 

dt  dx 


(i) 


(2) 


where  A  =  E~lA *,  B  =  E~1B‘.  If  A  has  only  real  eigenvalues  and  a 
complete  set  of  eigenvectors,  then  the  system  is  said  to  be  hyperbolic  (see, 
for  example,  Zauderer  [40]).  If  there  are  multiple  real  eigenvalues  and  less 
than  a  complete  set  of  eigenvectors,  then  the  system  is  of  (partial)  parabolic 


type.  If  all  of  the  eigenvalues  are  complex,  the  system  is  of  elliptic  type. 
Systems  with  complex  eigenvalues  are  not  causal.  Lyczkowski,  et  al.  [41], 
and  Sursock  [42]  provide  an  interesting  discussion  of  this  point  in  connection 
with  a  fluid  flow  problem.  The  underlying  problem  is  that  systems  with 
complex  eigenvalues  are  not  well-posed  as  initial  value  problems,  John  [43], 
Lax  [44].  We  will  not  consider  such  problems  any  further. 

If  E  is  singular,  (l)  can  give  rise  to  mixed  systems  of  all  types.  Some 
examples  can  be  found  in  Firedly  [45]  and  Lapidus  and  Pinder  [46].  Our 
interest  in  this  case  will  be  limited  to  purely  parabolic  systems  of  the  type 

(3) 

which  commonly  arise  in  engineering  problems. 


dw  _d2u;  ,  dw 

—  =  D——  +  A~~  +  Bw 
dt  dx 2  dx 


When  the  equations  of  motion  for  structural  elements  are  derived  from 
conservation  laws-in  particular,  from  variational  principles-the  resulting 
equations  are  typically  of  hyperbolic  type  (see,  for  example,  Crandall,  et  al. 
[38]).  However,  further  standard  assumptions  and  approximations  reduce 
the  equations  to  parabolic  systems  in  the  form  of  (3).  In  the  following 
paragraphs,  several  examples  will  be  given.  In  summary,  we  are  primarily 
interested  in  hyperbolic  systems  (l)  and  parabolic  systems  (3).  In  addition 
to  equations  (1)  or  (3),  there  are  associated  initial  and  boundary  conditions. 
For  equation  (1),  these  conditions  take  the  general  form 

initial  conditions  u>(x,  0)  =  f[x)  .  . 

boundary  conditions  Etu(0,  t)  +  Tw[L,t)  =  g(t)  ' 

where  dim(g)=dim(u;).  For  equation  (3),  they  take  the  general  form 


initial  conditions  tu(x,0)  =  /(x) 

boundary  conditions  Eju;(0,f)  +  +  Titi;(.L,t)  +  Tjj^iLyt)  =  g(t) 

(5) 


where  dim(y)=2dim(u;). 


It  is  well  known  that  the  coefficient  matrices  in  (4),  (5)  must  satisfy 
certain  constraints  if  the  problem  formulation  is  to  be  well-posed.  In  the 


hyperbolic  case  (equations  (1)  and  (4)),  these  constraints  essentially  require 
that  the  boundary  conditions  be  compatible  with  the  wave  directions.  Fur¬ 
ther  discussion  can  be  found  in  Russell  [30]  and  Agarwala  [47]. 


3.2  The  Timoshenko  Beam  Model 

We  will  show  how  some  conventional  beam  models  can  be  reduced  to  the 
standard  forms  described  in  the  preceding  paragraphs.  In  particular,  we 
will  begin  with  the  Timoshenko  model  and  then  consider  two  commonly 
used  approximations  which  can  be  derived  from  it,  the  Euler-Bernoulli 
model  and  the  “string”  model. 

Consider  the  beam  illustrated  in  Figure  1.  The  equations  of  motion  can 
be  derived  using  Lagrange’s  equations  and  in  the  absence  of  dissipation 
take  the  form 

"8  -  s  h' (S-*)]  W 

"8  ■  (£-*)] 

along  with  the  natural  boundary  conditions  for  a  =  0,  L 

displacement  or  shear  force 

r?(a,t)  =  qa(t)  kGA  _  <£(a,t))  =  fa{t)  ^ 

and 

rotation  or  moment 

4>(a,t)  =  El{^)=ra(t).  (8) 

The  two  equations  (6)  can  be  replaced  by  four  first-order  equations  by 
introducing  two  new  variables,  v[x,t)  and  7(x,t): 


dr) 

dv 

dt 

dx 

dv 

kG 

~di 

P 

w 


These  equations  clearly  represent  a  hyperbolic  system  and  the  natural 
boundary  c  onditions  become  for  a  =  0,L 


and 


displacement  or  shear  force 

l(o,<)  =  <?«,(!)  ■'(<>,<)  -  i/«(t), ■'.(<)  =  ^ 


(10) 


rotation  or  moment 

4>{a,t)  =  =  7 a(t),  7o(t)  =  ^ . 


(11) 


Note  that  the  boundary  conditions  applied  to  the  first  order  system 
(9)  require  the  time  integral  of  boundary  forces  or  moments  applied  to  the 
beam.  It  is  easy  to  confirm  that  the  transfer  functions  relating  forces  or 
moments  to  displacements  or  rotations  as  derived  from  either  equations  (6) 
or  (9)  are  indeed  identical  and  that  the  required  integration  is  essential. 


3.3  The  Bernoulli-Euler  Beam  Model 

The  Bernoulli-Euler  model  is  obtained  from  the  Timoshenko  model  with 
two  additional  assumptions: 

1.  rotational  inertia  is  neglected,  pi  — »  0 

2.  shear  deformation  is  neglected,  |3  —  <f>  — ►  0 


Assumption  (1)  reduces  the  second  of  equations  (6)  to 


Equation  (12)  and  assumption  (2)  are  now  used  to  reduce  the  first  equation 
of  (6)  to 

d 2  („Td2T)\ 
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Note  that  equation  (12)  along  with  assumption  (2)  leads  to  the  following 
expression  for  shear  force 


— (M --£("£)• 


Although  (14)  is  commonly  used  in  conjunction  with  the  Bernoulli-Euler 
model  (13),  it  should  only  be  used  with  caution.  Equation  (13)  is  valid 
only  in  the  limit  /  — +  0.  We  will  return  to  this  point  below. 

The  boundary  conditions  (7)  and  (8)  reduce  to  for  a  =  0,L 

displacement  or  shear  force 

•7(«,0  =1.(0.  =  /.(<),  (15) 


displacement  or  moment 

=  0  =*.(«),  ElS^a  =  r4t)  (16) 

Note  that  nonzero  shear  force  is  included  as  an  admissible  boundary  con¬ 
dition;  however,  the  remarks  following  equation  (14)  apply. 

Equation  (13)  can  be  reduced  to  first-order  form  by  introducing  a  new 
variable  7(1,  t) 

(17) 

dt  A  dx* '  [  ’ 

87  _  E  d2rt 
dt  ~  p  dx 2  ’ 

and  the  boundary  conditions  associated  with  (17)  are  for  a  =  0 ,  L: 


displacement  or 


shear  force 


iKO  =  „„((),  tfeu  =  v(,w  = 


and 


rotation  or 


^  =  *.w. 


moment 

iKO  =  ->.(<).  i.  =  ^ 


(19) 


Observe  that  (17)  is  a  parabolic  system  of  the  type  (3)  .  Equations  (17-  19) 
can  be  derived  directly  from  (13)  or  from  (9)  upon  invoking  assumptions 
(l)  and  (2).  We  should  also  note  that  a  corresponding  expression  for  shear 
force  obtains 


/ 


d 2/> 
dtdx 


(20) 


Because  shear  force  is  often  used  as  a  boundary  input  with  the  Bernoulli- 
Euler  model,  some  further  comment  is  warranted.  It  is  a  straightforward 
matter  to  compute  the  transfer  function  between  any  applied  boundary 
input  and  the  beam  response,  as  is  done,  for  example,  by  Kolousek  [48]. 
In  particular,  consider  the  deflection  response  at  x  =  0 ,H0[s),  to  a  force 
inputs  at  i  =  L,  Fl(s): 


H0() 


L 3  ^  cosh  A  cos  A  —  1  r  /  N 
A3(sinh  A  +  sin  A)  *L'S' 


(21) 


where 


A<(s)  =  - 


PAL\* 

El  ' 


The  transfer  function  in  (21)  has  a  branch  point  at  the  origin.  This  is 
true  of  any  transfer  function  associated  with  a  force  input.  On  the  other 
hand,  all  deflection  responses  to  torque  inputs  have  meromorphic  transfer 
functions.  Clearly,  there  is  an  issue  concerning  the  Bernoulli-Euler  model 
with  force  inputs  which  must  be  settled  before  it  can  be  used  with  any 
confidence. 


3.4  The  “String”  Model 


In  some  situations,  bending  deformation  may  be  negligible  with  respect  to 
shear  deformation,  that  is,  \4>\  C  \dr}/dx\.  In  this  case,  the  first  equation 


of  (6)  reduces  to 


with  boundary  conditions  for  a  =  0,  L 
displacement  or 


=  ya(t) 


shear  force 
KGA^=fa(t). 


(22) 


(23) 


This  simple  model  is  primarily  useful  for  illustrative  purposes.  Again  by 
introducing  the  new  variable  i/(z,t),  equation  (22)  can  be  replaced  by  two 
first  order  equations 


(24) 


dri 

dv 

dt 

dx 

dv 

KGdr i 

dt 

pdx 

which  is  to  be  solved  with  boundary  conditions  for  a  =  0,  L 

displacement  or  shear  force 

r)[a,t)  =  ria,  i/(a,t)  =  va(t),  va  =  kjp 


(25) 


Note  that  our  first  order  model  is  slightly  different  from  that  of  Burns  and 
Cliff  129). 


3.5  Distributed  Elements  with  Wave  Dynamics:  Fre¬ 
quency  Response  Calculations  for  State  Space 
Models 

We  will  be  concerned  with  the  computation  of  certain  irrational  transfer 
functions  and  a  resolvent  operator.  For  our  purposes,  the  res  -ent  can  be 
considered  as  an  integral  operator  with  kernel  called  the  Green's  function. 
In  this  section,  we  will  consider  explicitly  the  required  computations  for  the 
abstract  objects  discussed  previously.  To  do  so  we  will  focus  on  hyperbolic 
beam  models  for  distributed  elements.  Such  models  can  also  be  used  for 
elastic  dynamics  of  beams,  cables,  etc.  In  the  next  section,  we  will  discuss 


hybrid  system  models  consisting  of  interconnections  of  these  components 
with  rigid  bodies  and  other  lumped  parameter  models. 

Consider  a  class  of  hyperbolic  partial  differential  equations  in  one  space 
dimension  0  <  z  <  L,  arising  from  models  such  as  (??)  which  can  be 
written  (as  in  (2),  (4)) 

^^=A^^  +  Bx(i,2)  +  Cu(i.z)  (26) 

subject  to  boundary  conditions 

Ex(f ,  0)  +  Tx(t,L)  =  Df{t).  (27) 

Here,  x  is  an  n-vector  valued  state  x  E  Un(0,  L),  u  is  an  m-vector  val¬ 
ued  distributed  disturbance  u  (E  Mm(0,L),  A,B  are  real  n  x  n  matrices 
with  A  nonsingular  and  diagonalizable  [30],  and  E,  T  are  n  x  n  matrices. 
Controlability  questions  for  systems  of  this  type  are  considered  in  Russell 
[30].  We  remark  that  (28)  is  a  concrete  example  of  the  transform  of  the 
abstract  formula  (??).  Thus  it  is  clear  that  the  resolvent  for  the  operator 
A  :  >/"(0,L)  — ♦  ^"(0,  L)  defined  by  (26)  and  (27)  is  the  integral  operator 
Jo  Gr(s ,  2,  w)-dw.  After  taking  Laplace  transforms  in  the  temporal  variable 
t,  we  obtain  an  equation 

rl 

X{s,z)=  Gr{s,z,w)M[s,w)dw  +  J/flc(s,2).F(«)  (28) 

Jo 

where 

M(s,  w)  =  x(0,  w)  —  CU(s,  w), 

and  X,  U,  F  are  the  Laplace  transforms  of  x,  u,  /,  respectively.  The  function 
Gr(s,z,w)  is  the  Green’s  function  [49], [50]  for  (26),  (27)  and  Hbc{s,z)  is 
a  transfer  function  from  boundary  control  to  state.  Since  in  most  cases 
of  practical  interest  the  control  of  flexible  structures  will  be  effected  by 
actuators  whose  influence  functions  are  highly  localized,  we  have  formulated 
our  model  with  boundary  control  only. 


A  straightforward  calculation  leads  to  the  following  form  for  Hbc 

HBc{s,z)  =  *(s,z)[E  +  r*(s,L)]-1£>  (29) 


$(s,*)  =  el'1  (30) 

The  Green ’6  function  for  (26),  (27)  is  the  solution  to 

dGr{s^z,w)  =  A_,  [a/  _  Gr(5j2)U;)  +  U(z  -  w)  (3i) 

subject  to  the  boundary  conditions 

£Gr(s,0,w)  +  TGr(s,  L,w)  =  0  (32) 

where  £(•)  is  the  Dirac  delta  function  [49],  [50].  From  (31)  we  see  that  the 
solution  is  discontinuous  at  the  point  z  =  w.  After  some  computation,  we 
can  write 

f°r0<‘<"  (33) 

(  Gr  richt{s,z,w),  lot  w  <  z  <  L 


Gr  ieft{s,z,w)  =  Hbc[s,z)TQ{s,-w) 

Gr  right{s,z,w)  =  -Hbc{s,z)T9~1(s,L  -  w) 

3.6  Modeling  of  Hybrid  Systems 


In  most  applications,  models  for  the  dynamics  of  flexible  structures  will 
involve  interaction  between  various  elastic  and  rigid  elements.  In  the  par¬ 
ticular  case  of  flexible  structures  associated  with  large  space  structures, 
the  potential  topological  configurations  can  be  quite  complex.  Various  ele¬ 
ments  such  as  beams,  truss  structures,  cables,  membranes,  etc.,  may  have 
dominant  distributed  parameter  effects.  Typically  a  central  body  or  bodies 
represent  large  concentrations  of  mass  with  respect  to  the  overall  low  mass 
density  of  the  flexible  structure.  These  are  most  effectively  represented  by 
lumped  parameter  models  of  their  rigid  body  dynamics.  Additionally,  vari¬ 
ous  attitude  control  actuators  can  add  concentrated  inertia  elements  which 
can  be  effectively  modeled  as  lumped  Bystems.  Thus,  carefully  chosen  lin¬ 
ear,  hybrid  models  can  provide  an  effective  tool  for  analysis  of  dynamics 
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of  vibrations  and  their  effect  on  small  angle  motions  for  complex  space 
platforms.  In  this  section,  we  consider  the  structures  and  computations 
of  certain  resulting  transfer  functions  and  the  resolvent  operator  for  the 
composite  system  along  the  lines  of  Section  1. 

The  concept  of  a  mechanical  impedance  (terminology  borrowed  from 
electrical  network  theory)  has  been  used  in  structural  dynamic  modeling 
for  many  years  [48].  The  dynamic  stiffness  method  (application  to  space 
structure  modeling  is  reviewed  in  Piche  [51])  uses  this  notion  to  compute 
effective  transfer  function  models  for  interconnected  structures  [52].  Our 
approach  here  will  follow  along  similar  lines  except  that  we  will  focus  on 
computing  the  resolvent  operator  for  a  hybrid  structure  by  direct  manipu¬ 
lation  of  its  kernel;  viz,  a  Green’s  function,  xxx  A  hybrid  state  space  model 
is  constructed  in  Burns  and  Cliff  [29]  (where  considerations  are  given  for 
approximation  ad  computation  in  the  hybrid  state  space).  We  will  consider 
a  hybrid  state  space  as  consisting  of  a  direct  sum  of  spaces  X  =  X  <  ©  Xd 
where  X d  =  MNd  is  the  distributed  part  constructed  on  an  appropriate 
Hilbert  space  of  Nd-ve ctor  valued  functions  with  the  “energy”  inner  prod¬ 
uct  of  (??)  for  a  distributed  parameter  system  (DPS)  written  in  abstract 
form  (as  in  (??)) 

*)  +  Bu(0>  (36) 

with  x°d (z)  =  xrf(0,  z)  6  D{ A)  C  .  We  assume  that  A  generates  a  Co- 
semigroup  which  is  contractive  (so  that  (??)  is  well  defined).  We  ignore 
contribution  to  (??)  of  the  possible  distributed  control  force.  By  taking 
Laplace  transforms  in  (36),  we  can  write 

Xd[s,  z)  =  [si  -  A]~'BU[s)  (37) 

which  is  an  abstract  formulation  for  the  transform  of  (??). 

For  structural  control,  we  restrict  attention  to  the  hyperbolic  problem 
of  section  1.3: 

aZdgt'-  -  +  c„v(t,*)  (38) 

0  <  z  <  L,  subject  to  natural  boundary  conditions. 


Exrf(t,0)  +  Txd{t,L)  =  Dfd(t). 


mmm 
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For  this  problem,  equation  (37)  can  be  written 

Xd{s,z)=[  Gr{s,z,w)M{s,w)dw  +  HBc{s,q)Fd{s)  (40) 

Jo 

where 

M(s,  w)  =  Xj(w)  —  CV(s,  w),  (41) 

Ad,Cd,  Bd  are  matrices  defined  as  in  (2),  V  (s,  w)  is  the  Laplace  transform  o  f 
a  distributed  exogenous  disturbance,  and  Fj(s)  is  an  md- vector  of  inputs  to 
the  DPS.  Clearly,  (38)-(39)  can  represent  a  disjoint  collection  of  distributed 
elements  such  as  beams,  cables,  etc.  (Conceptually,  a  version  of  (40)  can 
also  be  written  for  higher  dimensional  spatial  domains,  but  we  feel  for  the 
current  presentation  that  the  required  complexity  of  notation  can  mask  the 
simplicity  of  the  underlying  concepts  (see  Butkovskiy  for  details  [50].) 

Similarly,  all  lumped  parameter  (LPS)  component  models  are  combined 
into  a  LPS  state  space  model  as 

xe(t)  =  Aexe(t)  +  x°t  =  xt( 0)  (42) 

with  x i  £  RNl  =  Xt  a  finite  dimensional  real  space.  By  taking  Laplace 
transforms  in  (42),  we  write  (analogous  to  (40)) 

Xt{s)  =  Rt{s)xi  +  H'{s)Ft{s ),  (43) 

where  Ri(s)  =  [s/^4  —  At ]-1  is  the  resolvent  for  the  (matrix)  operator  At 
and  Ht{s)  =  Rt{a)B. 

The  hybrid  state  space  X  =  Xt  ©  X  d  consists  of  elements 

xt{t) 
xd{t,z) 

which  are  N  =  Nd  +  iVrvalued  functions  of  z  £  [0,  L\,  t  >  0.  Finally, 
the  interconnection  of  component  systems  is  resolved  through  a  topological 
constraint  relation  consisting  of  m  =  md  +  mi  linear  equations; 


/(f)  +  K1xd(t,0)  +  K3xd(t,L)  +  K3xi(t)  =  Eu(t ) 


(45) 


where  u[t )  is  a  fc-vector  of  control  inputs  to  the  hybrid  system,  are 

m  x  Ni,  K3  is  m  x  Nt,E  is  m  x  k  real  matrices.  The  hybrid  modeling 
problem  is  to  find  an  equation  of  the  form  (40)  by  solving  (40),  (43)-(45) 
simultaneously  for  the  hybrid  state  x(t,z).  We  provide  the  resulting  model 
in  the  following  form: 

J\r(s,z)  =  f  Gr(s,z,w)M(s,w)dw  +  R(s,z)x°  +  H(s,z)U(s),  (46) 

Jo 

where  M(s,w)  is  given  in  (41).  The  resolvent  operator  for  the  hybrid 
system  is 

R(s;A)  —  f  Gr(s, z,w)  •  dw,  R(s,z )  (47) 

I/0 

where  R(s;  A)  :  X  — ►  Z?(A)  C  X ,  Gr  is  N  x  JV,*  and  R  is  N  x  Nt  are  matrix 
valued  functions  given  by 


i?(e  z\  —  Hi(s)Q1[s)  it  r  ( 

R{S'Z)  ~  l-2frc(«,*)QaW  J  3 ^(  )’ 


Gr(s,z,w)  = 


-Ife(s)Qi(s) 

Gr(s,z,w )  -  HBc{s,z)Q2{s) 


P{s,w) 


where 


<3M  =  |/«  +  «(.)|-‘  =  f9‘W 


j  >  v-y 

Q(s)  =  \K3Ht(S),  K^bcIs^  +  KiHbcIS'L)},  (51) 

P(s,w)  =  KiGr{s,0,  w)  +  K2Gr{s,L,w).  (52) 

Finally,  the  N  x  k  transfer  function  matrix  from  boundary  control  to  hybrid 
state  is 

*(«.«) -["f  *J^)]®(5)£-  (») 

The  derivation  of  (46)-(53)  is  straightforward  and  proceeds  as  follows. 
Substitute  (40),  (42)  into  (45)  and  solve  for  the  interconnecting  force  F(s). 
This  identifies  the  terms  Q(s),  P(s,w)  above.  Now  substitute  the  appro¬ 
priate  components  of  F(s)  into  (40),  (42)  and  use  the  hybrid  state  model 


4  Homogenization  of  regular  structures 


In  this  section  we  consider  the  problem  of  modeling  and  control  of  a  class 
of  lattice  structures,  e.g.,  trusses.  Their  large  size  and  repetitive  infras¬ 
tructure  require  special  techniques  for  structural  analysis  to  cope  with  the 
large  number  of  degrees  of  freedom.  Approximations  of  such  systems  by 
continua  provide  a  simple  means  for  comparing  structural  characteristics  of 
lattices  with  different  configurations,  and  they  are  effective  in  representing 
macroscopic  vibrational  modes  and  structural  response  due  to  temperature 
and  load  inputs.  Our  approach  to  the  construction  of  such  models  is  based 
on  a  technique  for  asymptotic  analysis  called  homogenization.  It  has  been 
widely  used  in  mathematical  physics  for  the  treatment  of  composite  systems 
like  porous  media  for  which  one  wishes  to  have  an  effective  approximating 
system  with  parameters  which  are  constant  across  the  structure.1 

Before  developing  the  general  features  of  the  method  and  applying  it  to 
the  treatment  of  lattice  structures,  we  shall  make  a  few  remarks  on  other 
work  on  continuum  models  which  has  appeared  in  the  recent  structural 
mechanics  literature. 

Noor,  et.  al.  [20]  use  an  energy  method  to  derive  a  continuum  approx¬ 
imation  for  trusses  with  triangular  cross  sections  in  which  the  modal  dis¬ 
placements  of  the  truss  are  related  to  a  linearly  varying  displacement  field 
for  an  equivalent  bar.  Plates  with  a  lattice  infrastructure  are  also  treated. 
In  Dean  and  Tauber  [9]  and  Renton  [19]  exact  analytical  expressions  for 
the  solutions  of  trusses  under  load  were  derived  using  finite  difference  cal¬ 
culus.  By  expressing  the  difference  operators  in  terms  of  Taylor’s  series 
Renton  [37]  was  able  to  derive  continuum  approximations  to  the  finite  dif¬ 
ference  equations  resulting  in  expressions  for  equivalent  plate  stiffnesses,  for 
example.  In  a  recent  paper  Renton  [37]  used  this  approach  to  give  equiv¬ 
alent  beam  properties  for  trusses,  which  complements  the  earlier  work  of 
Noor,  Anderson  and  Greene  [22],  and  Nayfeh  and  Hefzy  [21].  (See  also 
(Anderson  [22]).) 

^ee,  for  example,  the  papers  of  Larsen  [14],  Keller  [13],  and  the  reports  of  Babuska  [2] 
for  applications  and  discussions  of  design  techniques. 


37 


In  most  cases  a  continuum  model  is  associated  with  the  original  (lat¬ 
tice)  structure  by  averaging  the  parameters  of  the  lattice  over  some  natural 
volume  (e.g.,  of  a  “cell”  of  the  structure)  and  identifying  the  averaged  pa¬ 
rameter  value  (mass  density,  stress  tensor,  etc.)  with  the  corresponding 
distributed  parameter  in  the  continuum  model.  A  specific  form  for  the 
continuum  model  is  postulated  at  the  outset  of  the  analysis;  e.g.,  a  truss 
with  lattice  structure  will  be  approximated  by  a  beam,  with  the  beam  dy¬ 
namical  representation  assumed  in  advance.  While  this  approach  has  an 
appealing  directness  and  simplicity,  it  has  some  problems. 

First,  it  is  very  easy  to  construct  an  example  in  which  the  “approx¬ 
imate  model”  obtained  by  averaging  the  parameters  over  a  cell  is  not  a 
correct  approximation  to  the  system  behavior.  This  is  done  in  subsec¬ 
tion  4.1. 2  Second,  one  cannot  use  this  procedure  to  obtain  “corrections” 
to  the  approximation  based  on  higher  order  terms  in  an  expansion,  which 
may  sometimes  be  done  in  an  asymptotic  analysis.  These  terms  can  be  used 
to  describe  the  microscopic  behavior  (e.g.,  local  stresses)  in  the  structure. 
Third,  the  averaging  method  (averaging  the  parameters  over  space)  does 
not  apply  in  a  straightforward  way  to  systems  with  a  random  structure, 
since  the  appropriate  averaging  procedure  may  not  be  obvious.3  Fourth, 
the  method  cannot  be  naturally  imbedded  in  an  optimization  procedure; 
and  controls  and  state  estimates  based  on  the  averaged  model  may  not  be 
accurate  reflections  of  controls  and  state  estimates  derived  in  the  course  of  a 
unified  optimization  -  averaging  procedure.  In  particular,  the  method  does 
not  provide  a  systematic  way  of  estimating  the  degree  of  suboptimality  of 
controls  and  state  estimates  computed  from  the  idealized  model. 

In  this  work  we  use  a  totally  different  technique  called  homogenization 
from  the  mathematical  theory  of  asymptotic  analysis  to  approximate  the 
dynamics  of  structures  with  a  repeating  cellular  structure.  Homogenization 
produces  the  distributed  model  as  a  consequence  of  an  asymptotic  analysis 
carried  out  on  a  rescaled  version  of  the  physical  system  model. 

3See  the  numerical  experiments  in  (Bourgat  [7]). 

*  Homogenization  methods  do  apply  to  systems  with  a  randomly  heterogeneous  struc¬ 
ture,  see  (Papanicolaou  and  Varadhan  [16])  and  (Kunnemann  [39]).  We  shall  treat  such 
systems  in  a  subsequent  report. 


Unlike  the  averaging  method,  homogenization  cam  be  used  in  combina¬ 
tion  with  optimization  procedures;  and  it  cam  yield  systematic  estimates  for 
the  degree  of  suboptimality  of  controls  and  estimators  derived  from  ideal¬ 
ized  models.  While  our  results  are  stated  in  terms  of  simple  structures,  they 
demonstrate  the  feasibility  of  the  method;  and  they  suggest  its  potential  in 
the  analysis  of  structures  of  realistic  complexity. 

In  subsection  4.1  we  give  an  example  derived  from  (Bensoussan,  Lions, 
and  Papanicolaou  [5])  illustrating  some  of  the  subtleties  of  homogeniza¬ 
tion,  particularly  in  the  context  of  control  problems.  In  subsection  4.2  we 
derive  a  homogenized  representation  for  the  dynamics  of  a  lattice  struc¬ 
ture  undergoing  transverse  deflections.  We  show  that  the  behavior  of  the 
lattice  is  well  approximated  by  the  Timenshenko  beam  equation;  and  we 
show  that  this  equation  arises  naturally  as  the  limit  of  the  lattice  dynamics 
when  the  density  of  the  lattice  structure  goes  to  infinity  in  a  well  defined 
way.  The  problem  of  vibration  control  of  a  lattice  is  posed  and  discussed 
in  subsection  4.3.  In  subsection  4.4  we  derive  a  diffusion  approximation 
for  the  thermal  conductivity  of  a  one-dimensional  lattice  structure.  This 
property  is  useful  in  analyzing  new  materials  for  large  space  structures. 

Acknowledgements:  We  are  grateful  to  Professor  George  Papanico¬ 
laou  for  bringing  Kunnemann’s  paper  to  our  attention  and  to  Drs.  A.  Amos 
and  R.  Lindberg  for  their  comments  on  an  earlier  version  of  this  work. 

4.1  A  one-dimensional  example 

From  (Bensoussan,  Lions,  and  Papanicolaou  [5])  we  have  the  following  ex¬ 
ample: 


where  a*( i)  =f  a(x/c),  and  a(y)  is  periodic  in  y  with  period  Yo>  °(y)  >  &  > 
0.  It  is  simple  to  show  that 

IKIli.  =  /"  l«*(x)p  +  <  <  (2) 

and  so,  u‘  — ►  u  weakly  in  the  Hilbert  space  if1.4  Moreover, 

o-  -  M{a)  «  i  r°  a(y)iy  (3) 

Jo  ■'0 

and  it  is  natural  to  suppose  that  ue  — ♦  u  with  the  limit  defined  by 

-  -^[M(a)-^u(x)]  =  f(x)>x  e  ixo,Xi)  (4) 

u(x0)  =  ti(2l) 

This  is  untrue  in  general  (Bensoussan,  Lions,  and  Papanicolaou  [5],  pp. 
8-10).  The  correct  limit  is  given  by 

-  ^\a-^u(x)\  =  f(x),x  G  (xo,xi)  (5) 

u(x0)  =  u(xi) 

with 

5  «  [Af(i))->  (6) 

a 

In  general,  M(a)  >  a;  and  so,  the  error  in  identifying  the  limit,  (4)  versus 
(5),  is  fundamental. 

The  system  (4)  corresponds  to  averaging  the  parameter  a'(x)  over  a 
natural  cell;  a  procedure  similar  to  that  used  in  the  past  to  define  continuum 
models  for  lattice  structures.  As  (5)  shows,  the  actual  averaging  process 
can  be  more  subtle  than  one  might  expect,  even  for  simple  problems. 

4 Here  H 1  {«  €  Ls{*0,*i)  :  ||u||„,  <  oo) 


40 


£\ 


4.1.1  Homogenization  of  the  example 

To  see  how  (5)  arises,  we  can  use  the  method  of  multiple  scales  which  applies 
to  a  variety  of  perturbation  problems.  Suppose 

u‘(z)  =  u‘(x,  =  u0(x,  *)  +  £Ui(z,  *)  +  ...  (7) 

that  is,  we  suppose  that  u'  depends  on  the  “slow”  scale  x  and  the  “fast” 
scale  y  =f  x/e;  and  we  adopt  an  ansatz  which  reflects  this  dependence. 
Using  the  identity 

d  ,  .  x. ,  du  1  du  .  ,  x  ,  . 

—  u(x, With  y  = -.  (8) 

dx  e  ox  e  ay  e 

then  (l)  may  be  rewritten  as 

.  d  1  d  .  f  .  . .  d  1  <)  . .  r,  ,  ,  \ 

~  +  +  •••]>  =  f  (9) 

Simplifying  and  equating  coefficients  of  like  powers  of  e,  we  find  first  that 

‘  7^{a{y)^Uc]  =  °-  (10) 

The  assumptions  on  a(y)  imply 

«o(x,y)  =  u0(i)  (11) 

i.e. ,  no  y-dependence.  The  coefficients  of  e-1  satisfy 

.  d  .  .  d  ,  d  ...  d  .  d  ...  d  , 

{ai ;|a(!')siu“l  +  = 0  <12> 

or 

dy[  {V,dy  1J  dy  dx  {  ’ 

If  we  look  for  uj  in  the  form 


ui(z.y)  =  -x(y)-^  + 


(15) 


then  the  corrector  x(y)  must  satisfy 


“  T»ia{y)TyXMl  =  - 


da 

dy 


and  be  periodic.  That  is, 


a{y)zr  =  +  c 

dy 

(16) 

which  has  a  periodic  solution  (unique  up  to  an  additive  constant 

in  y)  if 

and  only  if 

1  c  , 

(17) 

—  1  +  r-rjdy  =  0 

Yo  Jo  a(y) 

which  implies 

c  =  @  a 

(18) 

We  obtain  an  equation  for  u0(x)  from  the  solvability  condition  for  u2(x,y). 
Equating  the  coefficients  of  c°  in  the  expansion,  we  have 


_^|o(!')|;“‘l "  = /w 
This  has  a  solution  u2(z,y),  periodic  in  y  if  and  only  if 


Jo  [«(v)  +  ^(«(y)x(y)]  -  a(y)~x(y)]dy}  •  -  (20) 

+/(*)  =  0 

where  we  have  used  (14).  The  integral  of  the  second  term  is  zero,  since  it  is 
the  integral  of  the  derivative  of  a  periodic  function  over  one  period.  Using 
(16)  and  (18),  (20)  reduces  to 


-  a 


d2u0 
dx 2 


+  f{x)  =  0 


(21) 


(plus  the  boundary  conditions)  which  is  (5)  (6). 


4.1.2  Control  and  homogenization  of  the  one  dimensional  sys¬ 
tem 


One  of  the  simplest  (stochastic)  control  problems  associated  with  the  pre¬ 
ceding  system  is  defined  by  the  Hamilton  -  Jacobi  -  Bellman  equation 


,x,d2u‘  1,-x.du'  .  ..1  ,  .  .  dul 

=  .el  2V  +  s(l'!')v7te  “ cu  1 


(22) 


x  €  0,u‘(x)  =  0  on  T  =  dO 

where  O  is  an  open  interval  in  R,  and  each  function  a(y),b{y ),  and  g(x,y) 
is  periodic  in  y  with  period  y0.  We  assume  that  a(y)  >  a  >  0  and  that 
c  >  0,  and  that  the  controls  v  take  values  in  R. 


This  Bellman  equation  corresponds  to  the  stochastic  control  problem 


u'(x)  =  inf  J*[v(-)] 
*(•) 


•/>(•)]  =  ex{  r  i(x',-,v)[r 

Jo  e 

dx'(f)  =  o(xf ,  — )dw(t )  +  -6(x‘,  —  )dt  +  G[x\  —,v)dt 
(  (  t  e 

x'(0)  =  xGO,l>0. 


(23) 


with  cr2(x,  y)  =f  a(y),  6(x,y)  =  b(y),  G(x,y,v)  =  g(x,y)v,  l{x,y,v)  =  \v2 , 
and  c(x,y,v)  =  c,  a  constant  in  (22).  Each  function  in  (23)  is  assumed 
to  be  periodic  in  y  with  period  one.  We  are  interested  in  the  behavior 
of  the  optimal  cost  and  control  law  for  (22)  in  the  limit  as  t  — *  0.  The 
stochastic  control  problem  (23)  was  treated  in  (Bensoussan,  Boccardo,  and 
Murat  [4]);  the  analysis  here  uses  different  arguments  which  emphasize  the 
computational  aspects  of  the  system. 


Evaluating  the  infimum  in  (22),  we  have  the  nonlinear  system 


a(“)u*x  +  7fe(7)u*  ~  cu*  ~  7)«)2  =  0 


t  t  it 

x  6  0,u'(x)|r  =  0. 


(24) 


The  analysis  of  the  control  problem  involves  homogenization  of  this  system. 


M  =  a{y)dvv  +  *>(y)dv 

with  its  formal  adjoint  defined  by 


A\  =  dv(a(y)dv-]  -  d„[(6(y)  -  a„(y))*]. 


The  problem 


A\m  =  0,y  — ►  m(y)  periodic  (27) 

m  >  0,  Jy  m(y)dy  -  1 

has  a  unique  solution  m(-)  on  V"  =  5°,  the  unit  circle,  with 

0  <  m  <  m(y)  <  m  <  oo.  (28) 

So  m(-)  is  a  density  on  Y.  We  cusume  that  b(-)  is  centered 

fY  ™(y)%)<*y  =  o.  (29) 

As  a  consequence  the  system 

Mx{y)  =  %)  (30) 

y  x(y)  periodic  ,J  x{v)dy  =  0 

has  a  well  defined  solution,  x(')  is  the  corrector  associated  with  the  problem. 
As  before  we  set  y  =  x/c  and  look  for  u'  in  the  form 


and  we  use 


u‘(x)  =  u‘(x,y)  =  u0(x,y)  +  euj(x,y)  +  ..., 


d,4>{x,y)  =  <fix{x,y)  +  ~4>v{x,y),y  =  x/f 
dzT<fi(x,y)  =  <f>xx(x,y)  +  2  f<£,„(x,y)  +  <t>vv(x,y). 
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Substituting  in  (24),  we  have 


a(y)[u0«*  +  2£Uoxy  +  ^t*oyy]  +  a(y)[tui*x  4-  2uiiy  +  -Uivv] 

H“a(y)  [e2u2*x  +  2 etl2zy  +  U2yy] 

+  ~&(y)[uOx  +  ~Uoy]  +  -6(y)[«uix  +  Uiy]  (33) 

+  *6(y)[£2u2l  +  eu2y]  -  c[u0  +  eui  +  £2u2] 

-^y2(x,y)[(u0l  +  £Ulx  +  £2u2l)  +  -(ti0y  +  ftily  +  f2«2y)P  =  0{  £2) 

The  last  term  is 

~  50a(**y)l(ittOy  +  2£U0lU0y  +  2u0xU!y  +  u2x)  (34) 

4-£(2uoiui1  +  2uxzUiy  +  2tijyU2y  -+-  2uoxu2v)]  +  0[  £2) 

Equating  coefficients  of  like  powers  of  £,  we  obtain 

(£-2Ka(y)uOyy  +  6(y)«0y  “  ^ {x,V)ulv  =  0}  (35) 

(f_1){a(y)uiyy  +  6(y)uiy  +  2a(y)u0*x  (36) 

+  i(y)u0*  -ffJ(l,y)«0*UOy  =  0} 

(f0){a(y)«2v»  +  &(y)«2y  +  2a(y)uity  +  fc(y)uu  (37) 

+a(y)u0lx  -  cu0  -  ^y2(x,y)u^  -  y2(x, y)u0xUix  =  0}. 

Choosing  u0(x,y)  =  u0(x),  which  must  be  justified,  satisfies  (35).  We  can 
then  solve  (36)  by  choosing 

«i(z»y)  =  -x(y)«o*(i)  +  «i(*).  (38) 

Equation  (37)  has  a  solution  for  u2(x,y)  if 

j  m(y){— 2o(y)xyuo«  -  6(y)xv«o*  +  a(y)uoxx 


(39) 


-euo  -  -0  ( x,y)(l  -  2xv)uoI}dy  =  °- 


This  gives  an  equation  for  u0(x) 


K..2 


qu0zz  -  cu0  -  ~Tu0x  =  0 


where 

q  d~  Jy  ~  2Xw (y)]  -  x(y)&(y)}<*y 

T  =f  j^m(y)g2(x,y)[l  -  2xv{y)}dy. 

Remark.  From  the  definition  of  Ai  and  the  corrector  x(y)  we  have 

fy  m{y)b{y)x{y)dy  =  JY\a(y)xvy  +  &(y)Xt,]x(y)™(y)<*y  (42) 

=  Jy  x(y)dvv[a(y)x(y)m(y)]dy  -  Jy  x(y)d|,[6(y)x(yMy)Hy 

Also,  using  (30), 

I  ™(y)%)x(y)<*y  =  j  x(y)a(y)m(y)xvv(y)dy  (43) 

j  y  jy 

-  Jy  x(y)fc(y)m(y)x»dy  +  2  Jy  x(y)xv^[a(y)"»(y)l«iy 

Adding  these  two  expressions,  we  have 

2  fYm{y)b(y)x(y)dy 

=  2  j  x(y)a(y)m(y)xvvdy  +  2  Jy  x(y)xv[om]vdy  (44) 

=  -2  /  dv\xam)xvdy  +  2  x(y)Xv\am\vdy  -  2  j  Xva(y)m(y)xvdy 

JY  JY  J  Y 

Thus,  q  may  be  rewritten  as 

9=  fYrn(y){a(y)\l-2xv  +  xl}}dy  (45) 

=  Jy  ^(y)a(y)[l  -  Xv?dy 


I 
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and  clearly  q  >  0. 

The  term  q  in  (41)  summarizes  the  effects  of  the  averaging  process 
on  the  uncontrolled  system.  The  homogenization  process  interacts  with 
the  control  system  through  the  term  T,  whose  form  would  be  difficult  to 
“guess”  from  simple  averaging  procedures. 

4.2  Continuum  Model  for  a  Simple  Structural  Me¬ 
chanical  System 

4.2.1  Problem  definition 

Consider  the  truss  shown  in  Figure  1  (undergoing  an  exaggerated  defor¬ 
mation).  We  shall  assume  that  the  truss  has  a  regular  (e.g.,  triangular) 
cross-section  and  no  “interlacing”  supports.  We  assume  that  the  displace¬ 
ments  of  the  system  are  “small”  in  the  sense  that  no  components  in  the 
system  buckle.  We  are  interested  in  describing  the  dynamical  behavior  of 
the  system  when  the  number  of  cells  (a  unit  between  two  (triangular)  cross 
sections)  is  large;  that  is,  in  the  limit  as 

€  =f  l/L  ->  0.  (46) 

We  shall  make  several  assumptions  to  simplify  the  analysis.  First,  we 
shall  assume  that  the  triangular  sections  are  essentially  rigid,  and  that  all 
mobility  of  the  system  derives  from  the  flexibility  of  the  members  con¬ 
necting  the  triangular  components.  Second,  we  shall  ignore  damping  and 
frictional  effects  in  the  system.  Third,  we  shall  confine  attention  to  small 
transverse  displacements  r ][t,x)  and  small  in  plane  rotations  4>{t,x)  as  in¬ 
dicated  in  Figure  1,  ignoring  longitudinal  and  out  of  plane  motions  and 
torsional  twisting.  Fourth,  we  shall  assume  that  the  mass  of  the  triangular 
cross  members  dominates  the  mass  of  the  interconnecting  links. 

Systems  of  this  type  have  been  considered  in  several  papers  including 
[20],  [21],  [22],  and  [37],  In  those  papers  a  continuum  beam  model  was 


Figure  1:  Deformed  truss  with  regular  cross-section. 


hypothesized  and  effective  values  for  the  continuum  system  parameters  were 
computed  by  averaging  the  associated  parameters  of  the  discrete  system. 
Our  approach  to  the  problem  is  based  on  homogenization  -  asymptotic 
analysis  and  is  quite  different. 

The  assumptions  simplify  the  problem  substantially,  by  suppressing  the 
geometric  structure  of  the  truss.  We  can  retain  this  structure  by  writing 
dynamical  equations  for  the  nodal  displacements  of  the  truss  members. 
For  triangular  cross  sections  nine  parameters  describe  the  displacements  of 
each  sectional  element.  The  analysis  which  follows  may  be  carried  over  to 
this  case,  but  the  algebraic  complexity  prevents  a  clear  presentation  of  the 
main  ideas.  As  suggested  in  (Noor  et  al.  [20])  one  should  use  a  symbolic 
manipulation  program  like  MACSYMA,  SMP,  Reduce,  etc.,  to  carry  out 
the  complete  details  of  the  calculations.  We  shall  take  up  this  problem 
on  another  occasion;  for  now  we  shall  treat  the  highly  simplified  problem 
which,  as  we  shall  see,  leads  to  the  Timoshenko  beam. 

We  shall  begin  by  reformulating  the  system  in  terms  of  a  discrete  el¬ 
ement  model  as  suggested  in  (Crandal  et  al.  [38]);  sec  Figure  2.  In  this 
model  we  follow  the  displacement  q,(f)  and  rotation  of  the  i,h  mass 

M.  The  bending  springs  (A:J)  tend  to  keep  the  system  straight  by  keeping 
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Figure  2:  A  lumped  parameter  model  of  the  simplified  truss  system. 


the  masses  parallel  and  the  shearing  springs  (k't)  tend  to  keep  the  masses 
perpendicular  to  the  connecting  links.  We  assume  small  displacements  and 
rotations  so  the  approximations 

sin^,(t)  «  &(t)  (47) 


tan  l\ni{t)H\  rii(t)/t 


are  valid. 


In  thi6  case  the  (approximate)  equations  of  motion  of  the  i1*  mass  are8 

<?4>*  .  l.<r,ffc+i(0  ~»7.(0, 


*,+i(0  - 


n,(t)  -  n.+i(0 


6 The  spring  constants  depend  on  *  since  they  represent  the  restorative  forces  of  flexed 
bars,  bent  by  different  amounts. 


where  we  have  normalized  M  —  1  and  defined 


si  =*  jlv-i  ~  Vi]  (50) 

and  similarly  for 

To  proceed,  we  shall  introduce  the  nondimensional  variable  e  =  tj L  and 
rewrite  the  system  (48)  (49)  as 

=  ijttvisw  -  «}  +  v+{x;v-^(i)}  (si) 

<f2n' 

where 

K\  =  k\l,Ki  =  k\l  (52) 

v<+rh  =  ^[r?,+1  -  ifc],  V'“ifc  =  i[r?<  -  rji _a]. 

Normalizing  t  —  1,  we  associate  a  position  x  G  [-1/2, 1/2]  with  each  mass; 
and  we  introduce  the  notation 

r]{t,Xi)  =  ri,(t),<i>(t,Xi)  =  <f>,{t).  (53) 

Having  normalized  t  =  1,  we  have  e  =  l  and  xl+1  =  x,  +  t  =  x,-  +  e.  Let 
Z  =  {x<}  be  the  set  of  all  points  in  the  system.  In  this  notation 

(V'+r7)(t,x)  =  *[r?(t,x  +  c)  -  r?(*,x)]  (54) 

(V‘~f?)(f,x)  =  j[r?(*,x)  -  T}[t,x  —  t)],x  G  Z 
and  the  system  is 

=  if.M{vV(i,xO  -  *•(!,*)> 

+  rV"{/f,(l.)V'-^(t,x,)}  (55) 

=  -V-{Jf.(x.)|VV(<,x,)  -  *•(! ,x,)]},x€  Z 


I 


The  scaling  of  (55)  may  be  interpreted  in  the  following  way:  Formally,  at 
least,  the  right  sides  of  both  terms  in  (55)  are  0(t-J).  This  implies  that  the 
time  variations  are  taking  place  in  the  “fast  time  scale”  r  =  t/e.  Also,  the 
spatial  variations  are  taking  place  in  the  “microscopic  scale”  x  which  varies 
in  e-increments  (e.g.,  x,+ j  =  x,  +  c).  Introducing  the  macroscopic  scale 
z  =  ex,  and  the  slow  time  scale  o  =  er,  we  may  rescale  (55)  and  observe 
its  dynamical  evolution  on  the  large  space-time  scale  on  which  macroscopic 
events  (e.g.,  “distributed  phenomena”)  take  place. 


Rewritten  in  this  spatial  scale,  the  system  becomes 


+  -ji-yjf, (*•)«•-*<(«,*•)} 


<F-'j,rV)  =  £«-{Jr.(V)[«‘+n(  <,*•)  -  «*•(<,  v)|} 


where 


6<±  =  fV«±  =  in  £ 


The  essential  mathematical  problem  is  to  analyze  the  solutions  4>e,ric  of 
(56)  (57)  in  the  limit  as  e  — ►  0. 


4.2.2  Mathematical  analysis 


To  proceed,  we  shall  generalize  the  problem  (56)  (57)  slightly  by  allowing 
K,  and  Kb  to  depend  on  z  as  well  as  z/c.  This  permits  the  restoring  forces 
in  the  model  system  to  depend  on  the  large  scale  shape  of  the  structure  as 
well  as  on  local  deformations.  We  use  the  method  of  multiple  scales;  that 
is,  we  introduce  y  =  x/e  and  look  for  solutions  of  (56) (57)  in  the  form 


r?‘(t)  =  T7  €(M,y), <£'(<)  =  4>'{t,z,V ), 


and  we  have 


K,  =  K,(z,y),Kb  =  Kb(z,y),y  =  - 


On  smooth  functions  rl){z,zf)  the  operators  6f±  satisfy 
(6t+xp)(z,y)  =  ip(z  +  e,y  +  1)  -  rp(z,y) 


=  ip{z,y  +  1)  -  ^(*,y)  +  ip{z  +  c,y  +  1)  -  xp{z,y  +  lj  (61) 

=  (S+0)(2,»)  +  y  +  1)  +  +  1)  +  0(<>) 

=  Wz*v)  -rp(z-€ty-  1) 

=  v>(*,y)  -  1>{z,V  -  1)  +  v{z,y  -  1)  -  Ip[z  -  c,y  -  1) 

=  (S"^)(s,y)  -  e|^(*,y  +  !)  +  ^f2f^'(2>s/  ”  J) +  °(f3)  (62) 

We  assume  that  <p(  and  rj f  may  be  represented  by 

4>f{t,z,y)  =  <t>o{t,z)  +  c<^i(f,«,y)  +  ...  (63) 


»?f(f,z,y)  =  t?o(*,*)  +  e*7i(M,y)  +  ... 

and  substituting  (63)  in  (56)  (57)  and  using  (60)  (61)  (62),  we  arrive  at  a 
sequence  of  equations  for  (<£o,  qo)>  [<f>i,  Vi),  ■  •  •  by  equating  the  coefficients 
of  like  powers  of  e. 

Starting  with  e_J, €_1,  e°, . . we  have 

^S+[rJf(z,y)S~0o(t.*)]  =  0  (64) 

which  is  trivially  true  from  (62)  (63).  The  same  term  involving  r/0(t,x) 
from  (57)  is  trivially  satisfied  by  the  assumption  (63).  Continuing 

j(S+{r/f6(z,y)5"<^i(f,z,y)}  (65) 

+K,{z,y){S+r}0{t,z)  -  <£o(M)}]  =  0 
which  may  be  solved  by  using  the  corrector  x*{ziV)  and  taking 

<t>i{t,z,y)  =  x4{z,y)<j>0{t,z)  (66) 
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at  ~  jkV.  /V  rfV  jA, 


S+{r.K:i(z,y)S  X*(*.y)}  =  K.{z, y)  (67) 

If  we  regard  z  as  a  parameter  in  (67),  then  there  exists  a  solution  x*>  unique 
up  to  an  additive  constant,  if  Ki,(z,  •),  K,  (z,  •)  are  periodic  in  y,  if  there  exist 
constants  A  said  B  so  that 

0  <  A  <  Kt(z,y)  <  B  <  oo  (68) 

and  if  the  average  of  K,(z,-)  is  zero 

1  fL /* 

jJ _l/2  K,(s,y)dy  =  0.  (69) 

Let  us  assume  that  (68)  (69)  hold,  and 

0  <  A  <  K,(z,y )  <  B  <  oo.  (70) 

Considering  (57),  the  0(e~l)  term  in  the  asymptotic  expansion  is 
,  ^[5~{/f4(z,y)(5+r/i(t,z,y)  -  0o(<, «))>]  =  0.  (71) 

Again  we  introduce  the  corrector  X^{z,y),  and  take  rji  in  the  form 

Vi{t,z,y)  =  Xn{z,y)4>  o{t,z)  (72) 

which  gives  the  equation  for  the  corrector 

5"{^.(2,y)[5+xn(2,y)  -  1]}  =  0  (73) 

or 

S-{K,{z,y)S+Xn{z,y)}  =  K,{z,y)  -  K.{z,y  -  1)  (74) 

By  hypothesis  the  right  side  in  (74)  is  periodic  in  y  and  has  zero  average 
(69).  Hence,  (74)  has  a  periodic  solution,  unique  to  an  additive  constant. 

Continuing,  the  O((0)  term  in  (56)  is 

S*  {rK6(z,y)S~ <t>2(t,z,y)}  +  K.{z,y)\S+r)i[t,z,y)  -  <t>i{t,z,y)} 


(75) 


+  K,[z,y)  —  rio (t.z)  +  S*{rKt(z,  y)  —  ^,(1.  2,  y)} 

Ql  Q  Q 

+  S+{rKb(z,y)  —  (f>o{t,z,y)}  +  —  {rA'6(z,y  +  1)} —<j>0{t ,  z) 

+§^irKdz,y  +  l)>^o(<.*)  -  =  o. 

This  should  be  regarded  as  an  equation  for  as  a  function  of  y  with  (t,z)  as 
parameters.  In  this  sense  the  solvability  condition  is  as  before,  the  average 
of  the  sum  of  all  terms  on  the  left  in  (75),  except  the  first,  should  be  zero. 
We  must  choose  (f>0  so  that  this  in  fact  occurs;  and  that  defines  the  limiting 
system. 

Using  the  correctors  (66)  (72),  we  must  have 
Average  ~  ^^[5+ (rift(z,y))  +  S+(rif6(z,y)x*(z,y))] 

~  ?£[£(rff.(*,V  +  1))]  -  (76) 

-M§^irKb{z^y  +  !))  +  S+{rKb{z,y)~xAz^y)) 

+K.(z,y)[S+Xr1{z,y)  -  x#(*,y))]}  =  0 

Defining  the  functions  EI(z),G(z)  by  the  associated  averages  in  (76),  the 
averaged  equation  is 

tst-r("wI 7>+cMfr-*M*  <77> 

which  is  the  angular  component  of  the  Timoshenko  beam  system  (Crandall 
et  al.  [38]  p.  348). 

Arguing  in  a  similar  fashion,  we  can  derive  the  equation  for  the  macro¬ 
scopic  approximation  displacement  of  the  lattice  system  in  terms  of  the 
“equivalent”  displacement  T]o(t,z)  in  the  Timoshenko  beam  system 

(78) 
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4.2.3  Summary 


We  have  shown  that  a  simplified  model  of  the  dynamics  of  the  truss  with 
rigid  cross  sectional  area  may  be  well  approximated  by  the  Timoshenko 
beam  model  in  the  limit  as  the  number  of  cells  (proportional  to  L/ i)  be¬ 
comes  large.  The  continuum  beam  model  emerges  naturally  in  the  analysis, 
as  a  consequence  of  the  periodicity  and  the  scaling. 

To  compute  the  approximate  continuum  model,  one  must  Bolve  (67)  and 
(74)  (numerically)  for  the  correctors  and  then  compute  the  parameters  in 
(77)  (78)  by  numerically  averaging  the  quantities  in  (76)  (and  its  analog 
for  (57))  which  involve  the  correctors  and  the  data  of  the  problem. 


4.3  Homogenization  and  Stabilizing  Control  of  Lat¬ 
tice  Structures 

In  this  subsection  we  show  that  the  process  of  deriving  effective  “contin¬ 
uum”  approximations  to  complex  systems  may  be  developed  in  the  context 
of  optimal  control  designs  for  those  systems.  This  procedure  is  more  effec¬ 
tive  than  the  procedure  of  first  deriving  homogeneous  -  continuum  approx¬ 
imations  for  the  structure,  designing  a  control  algorithm  for  the  idealized 
model,  and  then  adapting  the  algorithm  to  the  physical  model.  In  fact, 
separation  of  optimization  and  asymptotic  analysis  can  lead  to  incorrect 
algorithms  or  ineffective  approximations,  particularly  in  control  problems 
where  nonlinear  analysis  (e.g.,  of  the  Bellman  dynamic  programming  equa¬ 
tion)  is  required. 

We  shall  apply  the  combined  homogenization  -  optimization  procedure 
described  in  subsection  4.1  (based  on  (Bensoussan,  Boccardo,  and  Murat 
[4]))  to  the  problem  of  controlling  the  dynamics  of  lattice  structures  like  the 
truss  structure  analyzed  in  the  previous  subsection.  We  shall  only  formulate 
a  prototype  problem  of  this  type  and  discuss  its  essential  features. 

Consider  the  model  for  the  lattice  structure  analyzed  in  subsection  4.3 


Figure  3:  Truss  with  transverse  actuator  forces. 


with  control  actuators  added.  The  truss  shown  in  Figure  1  iB  again  con¬ 
strained  to  move  in  the  plane  and  torsional  motion  is  excluded  to  simplify 
the  model  and  confine  attention  to  the  basic  ideas.  Now,  however,  we  in¬ 
clude  a  finite  number  of  actuators  acting  to  cause  transverse  motions.  The 
truss  with  actuator  forces  indicated  by  arrows  is  Bhown  in  Figure  3.  The 
corresponding  discrete  element  model  is  shown  in  Figure  4. 

Suppose  that  the  physical  actuators  act  along  the  local  normal  to  the 
truss  midline  as  shown  in  the  figures,  and  that  the  forces  are  small  so  that 
linear  approximations  to  transcendental  functions  (e.g.,  sin^  «  ^>i,  etc.) 
are  valid.  Then  the  controlled  equations  of  motion  of  the  discrete  element 
system  are  (recall  equation  (51)) 

rd-^T  =  ^:{V'+,;(0  -  <«(()}  (79) 

-3-  -  -  *■(!)]>  +  £*(•',•>,(<) 

/=! 

where  the  notation  in  (52)  has  been  used, 
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Figure  4:  Discrete  element  model  of  the  controlled  truss. 

and  t;,  j  =  1, . . .  ,m  are  the  locations  of  the  actuators.  Hence,  if  =  0 

for  all  j  =  1, . . .  ,m  there  is  no  actuator  located  at  the  ith  point  which  corre¬ 
sponds  to  the  physical  point  x  €  [0,  L\.  The  number  m  of  actuators  is  given 
at  the  outset  and  does  not,  of  course,  vary  with  the  scaling.  The  control 
problem  is  to  select  the  actuator  forces  as  functions  of  the  displacements 
and  velocities  of  components  of  the  structure  to  damp  out  motions  of  the 
structure.  Measurements  would  typically  be  available  from  a  finite  number 
of  sensors  located  along  the  structure.  We  shall  not  elaborate  on  this  com¬ 
ponent  of  the  model,  and  shall  instead  assume  that  the  entire  state  can  be 
measured.  To  achieve  the  stabilization,  we  shall  associate  a  cost  functional 
with  the  system  (79).  Let 

u(0  ~  [ui(0* •  •  •  *um(0]T  (81) 

be  the  vector  of  control  forces,  and 

•'XH  =  /"£<«■.  W(<)]’+Mo:(‘)]’ 

^  ~  a  t j.* /#\1  J  I o*>\ 


+  D«(Mi)u*(0>e-** 

i= i 

where  (a,, 6.)  and  (a,,/?,)  are  non-negative  weights.  Formally,  the  control 
problem  is  to  select  =  l,...,JV,j  =  1  to  achieve 

inf  •HOI  (83) 

subject  to  (79)  (80)  and  the  appropriate  boundary  conditions.  The  case 
"7  — *  0  corresponds  to  stabilization  by  feedback. 

The  analysis  of  this  control  problem  is  based  on  the  scaling  used  in 
subsection  4.3,  equations  (51)  -  (58).  Let  r  =  t/t  be  the  fast  time  scale, 
then 

J1l«(-)]  =  /  «X!W&(t)]2  +  Mft(r)]2 

Jo  <=i 

+  a,cJ[^,(r)]2  +  ft«2(^'(r)]2  (84) 

+  £6(‘>*>)u2(r)}c~t7rrfr 

;= i 

with  </>'(r)  =  <£'(cr),  etc. 

Let  (<f>,  rj,  17)  be  the  state  vector  of  the  system  (79)  with  <£  =  [^1, . . . ,  <f>n 
and  similarly  for  the  other  terms.  Let  V  =  Vt'1(<l>,^>,r],r])  be  the  optimal 
value  function  for  the  problem  (79)  (84).  Then  the  Bellman  equation  asso¬ 
ciated  with  (79)  (84)  is 

< 

t=i 

t=l  r  r 

+  (f;{-V'-[ir;(v-,,-^,)))v,,  (85) 


££!*(•'. -.Kv*  +  <(>,i,)uj|> 


+£  £3(a,4>,2  +  M,2  +  e2(a,^2  +  /*.»? ,2)]  -  CTfV  =  0. 

•=i 

Remarks: 

1.  Note  that  the  minimization  in  (85)  is  well  defined  if  the  admissible 
range  of  the  control  forces  is  convex  since  the  performance  measure 
has  been  assumed  to  be  quadratic  in  the  control  variables  £(i,t;)u;. 

2.  Since  we  have  not  included  the  effects  of  noise  in  the  model,  the 
state  equations  are  deterministic  and  the  Bellman  equation  (85)  is  a 
first  order  system.  To  “regularize”  the  analysis,  at  least  along  the 
lines  followed  in  conventional  homogenization  analysis,  it  is  useful  to 
include  the  effects  of  noise  in  the  model  and  exploit  the  resulting 
coercivity  properties  in  the  asymptotic  analysis. 

3.  If  we  introduce  the  macroscopic  spatial  scale  z  =  ex,  the  mesh  {x,}, 
and  the  variables 

<f>{t,z,)  =  <j>\{t),4>{t,Zi)  =  #(t),  etc.  (86) 

then  the  sums  may  be  regarded  as  Riemann  approximations  to  inte¬ 
grals  over  the  macroscopic  spatial  scale  z.  The  asymptotic  analysis  of 
(85)  with  this  interpretation  defines  the  mathematical  problem  con¬ 
stituting  simultaneous  homogenization  -  optimization  for  this  case. 

4.4  Effective  conductivity  of  a  periodic  lattice 

In  this  subsection  we  consider  a  version  of  a  heat  conduction  problem 
treated  by  Kunnemann.  Simple  expressions  for  thermal  properties  of  com¬ 
posite  materials,  have  been  derived  in  the  past  using  homogenization  tech¬ 
niques.  The  derivation  of  effective  conductivities  for  discrete  structures 
is  useful  for  assessing  the  thermal  response  of  such  structures  in  variable 
environmental  conditions. 
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4.4.1  Problem  definition 


Let  Z  —  {0,  ±1,  ±2, . . .}  and  Zd  =  Z  x  . . .  x  Z  ( d  times)  be  a  d-dimensional 
lattice.  Let  (  >  0  be  a  number  small  relative  to  1.  We  want  to  describe 
the  effective  conduction  of  thermal  energy  on  the  e-spaced  lattice  eZd.  Let 
e,  =  (0,0, ...  0, 1,0, .. .  0)r  with  1  in  the  ith  position,  i  =  1,2,  . . .  ,d.  If  x  is 
a  point  in  eZd,  then  z  ±  ee,,  1  <  t  <  d,  are  the  nearest  neighbors  of  x.  Let 
a±(x),x  £  Zd,  1  <  i  <  d,  be  the  two  functions  defined  on  the  lattice,  and 
assume 


a,(z)  d=  a+(z)  =  dt-(x  -t  e,),z  £  Zd,  1  <i  <d 
0  <  A  <  a,-(z)  <  B  <  oo,  Vz  £  Zd,  1  <  t  <  d 
a,(z)  is  periodic  with  period  t  >  1 
in  each  direction,  1  <  i  <  d  * 

Next  let 

ai±(x)  =  €  eZd,  1  <i<d. 


(87) 

(88) 
(89) 


(90) 


Equation  (88)  means  that  the  conduction  process  is  reversible  and  that 
the  conductivity  Oi(z)  is  a  “bond  conductivity,”  i.e.,  independent  of  the 
direction  in  which  the  bond  (z,z  +  e,)  is  used  by  the  process.  Equation 
(90)  means  that  the  configuration  of  bond  conductivities  a'±(-)  on  eZd  is 
simply  a,±(-)  on  tZd  “viewed  from  a  distance.”  Assumption  (89)  imposes  a 
regularity  condition  on  the  physics  of  the  conduction  process.  An  assump¬ 
tion  like  this  is  essential  for  existence  of  a  limit  as  e  — ♦  0.  In  one  dimension 
the  situation  is  illustrated  in  Figure  5  and  Figure  6.  A  system  similar  to 
this  with  random  bond  conductivities  was  treated  by  Kunnemann  [39]  by 
imposing  some  ergodicity  properties  on  the  bond  conductivities. 

One  can  associate  with  this  system  a  random  (jump)  process 

{*f(t,z),t  >  0,x  £  iZd) 


®The  period  may  be  different  in  different  directions. 


on  the  e-spaced  lattice.7  In  effect,  as  t  — *  0  ,  {X‘}  converges  to  a  Brownian 
motion  on  the  lattice;  and  the  main  result  of  the  analysis  is  an  expression  for 
the  diffusion  matrix  Q  =  [gtJ;  i,j  =  1, 2, . . . ,  d]  of  this  process.  This  matrix 
describes  the  macroscopic  diffusion  of  thermal  energy  in  the  system.  It  is 
the  effective  conductivity. 

We  shall  carry  out  the  asymptotic  analysis  of  this  system  in  the  limit 
as  e  — +  0  using  homogenization.  Let 

(Vpu)(x)  d=  j[u(x  -  ce.)  -  u(x)j  (91) 

(V'+u)(x)  d=  i[u(x  +  ee.)  -  u(x)] 
x  G  f  Zd,  1  <  i  <  d, 

for  any  u  square  summable  on  eZd  or  square  integrable  on  Rd  with  e,  the 
ith  natural  basis  vector  in  Rd.  Then 

ol  «=i  c 

=  LV(t,x) 

is  the  diffusion  equation  on  the  t-spaced  lattice  with  density  u'(f,x)  and 
conductivity  Oj(x/f).  We  are  interested  in  an  effective  parameter  represen¬ 
tation  of  the  thermal  conduction  process  as  e  — *  0. 

Remark:  Although  probabilistic  methods  are  not  required  in  the  anal¬ 
ysis,  the  associated  probabilistic  framework  has  a  great  deal  of  intuitive 
appeal.  The  operator  U  may  be  identified  as  the  infinitesimal  generator  of 
a  pure  jump  process  X‘(s)  in  the  “slow”  time  scale  s  =  e2t;  (Breiman  [8]). 
Moreover,  Ll  is  selfadjoint  on  tZd  with  the  inner  product 

(/><?)  =  J2  /(*)»(*)•  (94) 

*ez* 


7 Definition  of  this  process  is  not  necessary  for  the  analysis,  but  it  bolsters  the  intuition. 


Hence,  the  backward  and  forward  equations  for  the  process  A'*  (a)  are,  re¬ 
spectively, 

dp  (iMl1)  _  |r<  1/  ,  .1  ir.r\ 


dp‘(y,<|x) 


=  (I‘p‘(y,i|-)!(x) 

=  [IV(-,(|x)]M 


So  the  process  is  “symmetric”  in  the  sense  of  Markov  processes  (Breiman 


The  asymptotic  analysis  of  (93),  when  interpreted  in  this  context,  means 
that  as  the  bond  lattice  is  contracted  by  c  and  time  is  sped  up  by  f-2, 
the  jump  process  {A”‘(s)}  approaches  a  diffusion  process  with  diffusion 
matrix  Q.  In  other  words,  on  the  microscopic  scale  thermal  energy  is 
transmitted  through  the  lattice  by  a  jump  process;  but  when  viewed  on  a 
macroscopic  scale  the  energy  appears  to  diffuse  throughout  the  lattice.  The 
microscopic  physics  are  described  in  (Kirkpatrick  [1 1])  and  (Kittel  [12]). 
The  approximation  developed  belo  w  for  a  periodic  lattice  is  similar  to 
the  one  developed  by  Kunnemann  for  a  random  lattice.  This  similarity 
demonstrates  the  robustness  of  the  method,  and  the  limited  dependence  of 
the  macroscopic  properties  of  the  medium  on  the  details  of  the  microscopic 
variations  of  the  structure. 

Because  the  basic  problem  (93)  is  “parabolic,”  we  can  introduce  the 
probabilistic  mechanism  and  make  use  of  it  in  the  analysis.  In  the  “hyper¬ 
bolic,”  structural  mechanical  problems  we  treated  before  this  device  was 
not  readily  available. 


4.4.2  Asymptotic  analysis-homogenization 


The  essential  mathematical  step  is  to  show  strong  convergence  of  the  semi¬ 
group  of  L€,  say 


r  <*«f  r«i  rrh\  d«f  Lt 


and  to  identify  the  limiting  operator 


d 2 

3  dxidij 

This  is  accomplished  by  proving  convergence  of  the  resolvents 
,  [  — 2/  +  a]-1 — *\—L  +  a]-1  as  e  — ♦  0 
That  is,  if  /  is  a  given  function  and 

«*(•)  =  [-!*  +  a]'1/ 

«(.)«[-X +  *]->/ 

then  u‘  -*  u  in  an  appropriate  sense. 


(97) 


(98) 

(99) 


The  method  of  multiple  scales  will  be  used  to  compute  the  limit.  Be¬ 
cause  the  conductivities  at(x)  in  (93)  do  not  depend  on  time,  we  may  work 
directly  with  Ll  rather  than  the  parabolic  PDE  (93)  (cf.  (Bensoussan,  Li¬ 
ons,  and  Papanicolaou  [5])  Remark  1.6,  p.  242).  The  method  of  multiple 
scales  is  convenient  because  it  is  a  systematic  way  of  arriving  at  the  “right 
answers”  -  something  which  is  not  always  simple  in  this  analysis. 

Bearing  in  mind  (99),  we  consider 

(LV)(x)  =  /(x)  (100) 

with  u'(x)  in  the  form 

u‘(x)  =  u0{x,  *)  +  cuj(x,  *)  +  e2u2(x,  *)  +  . . .  (101) 

with  the  functions  u,(x,  y)  periodic  in  y  €  f.Zd  for  every  j  =  0, 1,. . ..  (As  it 
turns  out  the  boundary  conditions  are  somewhat  irrelevant  to  the  construc¬ 
tion  of  “right  answers.”)  To  present  the  computations  in  a  simple  form,  it  is 
convenient  to  introduce  y  =  x/c,  to  treat  x  and  y  as  independent  variables, 
and  to  replace  y  by  x/c  at  the  end. 

Recall  the  operators  V'*  from  (93).  Applied  to  a  smooth  function 
u  =  u(x,x/c),  we  have 

(V'"u)(x,y)  =  -[tt(x  -  ce„y  -  e<)  -  u(x,y)] 


(102) 


=  -[u(z,y  -  d)  -  u(x,y)]  +  -[u(x  -  e.)  -  u(z,y  -  e,)j 

=  ^(Vr«)(z,y)-  |^(x,y-e,)  +  £^^(x,y-e,)-l-0(£3) 
where  on  functions  (f>  —  <f>(y) 

(v,~ 4>){y)  =  4{u  ~  c«)  -  Hy)  (103) 

Defining 

(v,+^)(y)  =  4>[y  +  «<)  -  4>{y)  (104) 

we  also  have 

(V'+u)(x,y)  =  j(V,+u)(x,y)  +  ^-(x,y  +  e<)  (105) 

1  d2u .  .  . 

+e2a®^*,y  +  e^  +  °^’ 


Now  we  substitute  (101)  into  (100)  and  use  the  rules  (103)  (104).  Equat¬ 
ing  coefficients  of  like  powers  of  e,  this  leads  to  a  sequence  of  equations  for 
u0,ui, . . ..  Specifically,  (using  the  summation  convention) 

(Z/u')(x,y)  =  -Vr[o,(y)V'+«'] 

=  ^vrb(y)v,+uo(*>y)]  ~  vrMy)|—(*>y  +  «0] 


~\€Vi  [a»(y)^(l.y  +  e01  +  c>(e) 

-jvr[o<(y)v<+tt1(x,y)] 

-fV‘-[ai(y)^i(x,y  +  e,)]  +  O(f) 

-v.r[«.-(y)v<+«2(*,»)]  +  o(0  =  /(x) 
That  is,  labeling  each  term  by  its  order  in  £ 

(£-’)  V-[a,(y)V  +  u0}  =  0 


(106) 


(107) 


imm 

CfR 

(«  )  «v,‘  [“i(y)^7(i.y  +  «.)1  +  v.  k(y)vi  «i(i»y)l  =  0  (iog) 

and  (recall  cV'*  is  0(1)  in  c) 

(t°)  ^£V‘-(o,(y)^^(i,y  +  e,)] -<V'~(o,(y)^(i,y  +  e,)]  (109) 

-V,"[G1(y)Vl+u2(i,y)]  =  /( x) 

From  (107)  we  have 

o,(y  -  e,)lu0(i,  y)  -  U0(i,y  -  e*)]  (110) 

-a.(y)[«o(x,y  +  ej)  -  u0(x,y)]  =  0 

If  we  take  u0(z,y)  =  u0(x),  this  is  trivially  true;  and  (108)  simplifies  to 


eV'  +  V<  My)V,tui(*>y)]  =  °- 

At  this  point  we  introduce  “correctors,”  That  is,  we  assume 

*M*.y)  =  Y  x*(y)§^  +  M*) 

*=1 


(111) 


(112) 


with  x*(‘)  the  correctors.  Using  this  in  (111),  we  have  (again  using  the 
summation  convention) 

^r[a.(y)v,+Xfc(y)]^  +  [a*(y-«*)-ot(y)]|^  =o  (113) 

If  we  take  x*(y)  48  the  solution  of 


V,  My)V,+  Xk(y)|  +  [a*(y  -  Ci)  -  a*(y)]  =  0 


(114) 


(we  have  to  verify  the  well-posedness  of  (114)),  then  (113)  is  satisfied.  (The 
term  uj(z)  is  determined  (formally)  from  the  O(e)  term  in  the  system  (101) 
(106).) 


Regarding  the  well-posedness  of  (114),  note  that 

^■K(y)v(**(»)|  =  V’(v) 


(115) 


has  a  periodic  solution  on  tZ  which  is  unique  up  to  an  additive  constant  iff 
the  average  of  the  function  V>(y)  over  a  period  (c£)  is  zero;  i.e., 


V’=  +  M  =0,n  =  l,2,...,<f.  (116) 

L  *=i 

This  condition  clearly  holds  in  (114),  and  so,  Xt(y)  is  well  defined  (up  to 
an  additive  constant). 


We  shall  determine  the  equation  for  uo(x)  by  using  (112)  (114)  in  (109). 
Using  the  Kronecker  delta  function  6ik,  we  have 


l«.(y)**] 


d2u0 

dztdxk 


K(y)xk(y +  «.-)) 


d2u0 

dxtdxk 


/(*) 


=  {^V-[o,(y)^]  -  V-[a,(y)V+X*]}^|^  (117) 

-VrKtyJXkty)]^^  -  Vf  [a,(y)V+u2]  =  /(*) 

The  term  in  braces  is  zero  from  (114).  To  obtain  the  solvability  condition 
(116)  for  u2,  we  introduce  the  average 


=f  symmetric  part  {-V~{a,(y)x*(y)]} 

Then  solvability  of  the  equation  for  u2  gives  the  equation 

d2u0  t,  N 
=  /w' 


(118) 


(119) 


And  this  is  the  diffusion  equation  which  defines  the  limiting  behavior  of  the 
system  (100)  in  the  macroscopic  i-scale  in  the  limit  as  c  — ♦  0. 


We  can  justify  the  asymptotic  analysis  by  using  energy  estimates  or 
probabilistic  methods  as  in  (BenBoussan,  Lions,  and  Papanicolaou  [5]).  (See 
also  Kunnemann  [39]).)  We  shall  omit  this  analysis  here. 
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4.4.3  Summary 


Returning  to  the  original  problem  (93)  for  the  evolution  of  thermal  energy 
on  a  microscopic  scale,  we  have  shown  that  the  thermal  density  u'(t,x)  — ♦ 
u0(*,:r)  as  f  — ►  0  (in  an  appropriate  norm)  where 


du0  _  1  y-v  du0 
dt  2  |/“1  dx,dij 

with 

=  -j  XHV.+My)x*(y)]  +  k(y)x,(y)]} 

L  k=l 

with  the  correctors  Xk,k  =  1,2 given  by 

Y,  v.rMv)v<+x<(v)]  =  -la*(y  - e*)  -  ak{y)  1 

•=i 


(120) 


(121) 


(122) 


k  =  1,2,. ..  ,d 

To  compute  the  limiting  “homogenized”  model  (120),  one  must  solve  the 
system  (122)  (numerically)  and  then  evaluate  the  average  (121). 


The  fact  that  the  original  problem  (93)  is  “parabolic”  (i.e.,  it  describes 
a  jump  random  process),  enables  us  to  exploit  the  associated  probabilistic 
structure  to  anticipate  and  structure  the  analysis.  In  this  way  we  can 
anticipate  that  the  limit  problem  will  involve  a  diffusion  process.  T"  fact, 
the  arguments  used  are  entirely  analytical  and  the  limiting  diffusion  (120) 
is  constructed  in  a  systematic  way.  It  is  not  postulated.8 


*Prob&bili«tic  arguments  can  be  used  (Bensoussan,  Lions,  and  Papanicolaou  [6],  Chap¬ 
ter  3);  and  they  have  some  advantages. 
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